In [1]:
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import mapclassify
In [ ]:
 
In [2]:
def wf(x):
    return "/Users/andywang/desktop/Lab3/" + x
shp = gpd.read_file(wf("pc_MA_EE508.shp"))
d = pd.read_csv(wf("pc_MA_EE508.csv"))
shp_ma_outline = gpd.read_file(wf("MA_outline.shp"))
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/3228652068.py:4: DtypeWarning: Columns (20) have mixed types. Specify dtype option on import or set low_memory=False.
  d = pd.read_csv(wf("pc_MA_EE508.csv"))
In [3]:
shp = shp.set_index('pid')
d = d.set_index('pid')
In [4]:
shp_c = shp.copy()
shp_c['geometry'] = shp_c['geometry'].centroid 
In [5]:
shp_c['ms'] = d['ha'] / 10
shp_c.join(d['travel']).plot(column='travel', markersize='ms', alpha=0.5, scheme = 'quantiles')
Out[5]:
<Axes: >
In [6]:
# Compute the change in percentage of parcel-level forest cover between 1990 and 2010
d['p_f_change'] = d['p_f_2010'] - d['p_f_1990']

# Compute the log of the per-hectare sale price of the last sale of the property
d['ls_price_per_ha_log'] = np.log(d['ls_price'] / d['ha'])

# Convert the 'ls_date' column to a datetime format and extract the year
d['ls_year'] = pd.to_datetime(d['ls_date']).dt.year

# Create a new variable 'p_e_ls' to represent the percentage of easement at the time of the last sale
# If 'e_year' is less than or equal to 'ls_year', use 'p_e'; otherwise, set it to 0
d['p_e_ls'] = np.where(d['e_year'] <= d['ls_year'], d['p_e'], 0)

# Normalize the building area to the parcel size
d['bld_area_per_ha'] = d['bld_area'] / d['ha']
In [7]:
# Define the list of predictor variables
XCOLS_LS = ['slope', 'travel', 'p_wet', 'coast_2500', 'river_frontage', 
            'lake_frontage', 'pop_dens_tract_2012-2016', 
            'hh_inc_avg_tract_2012-2016',
            'p_f_2000', 'p_prot_2000_500', 
            'p_d_2000_200', 'p_d_2000_1000', 'p_d_2000_5000',
            'p_dev_open', 'p_dev_low', 'p_dev_medium', 'p_dev_high', 
            'p_crops', 'p_pasture','p_grassland', 'bld_area_per_ha', 
            'ls_year']

# Create the X and y DataFrames/Series for last sale price prediction
X_ls = d[d['ls_year'] >= 2000][XCOLS_LS]  # Select data from 2000 and onward
y_ls = d[d['ls_year'] >= 2000]['ls_price_per_ha_log']  # Sale price per hectare with log transformation
In [8]:
# Define the list of predictor variables
XCOLS_FC = ['slope', 'travel', 'p_wet', 'coast_2500', 'river_frontage', 
            'lake_frontage',  'pop_dens_tract_2012-2016', 
            'hh_inc_avg_tract_2012-2016', 
            'p_f_1990', 'p_prot_1990_500', 
            'p_d_1990_200', 'p_d_1990_1000', 'p_d_1990_5000']

# Create the X and y DataFrames/Series for forest change probability estimation
X_fc = d[XCOLS_FC]  # All data without any time filter
y_fc = d['p_f_change']  # Change in percentage of forest cover from 1990 to 2010
In [9]:
# Define the selection criteria for X_ls and y_ls
i_ls = d[XCOLS_LS].isnull().sum(1).eq(0)
i_ls &= d['ls_year'].ge(2000)
i_ls &= d['ls_year'].le(2019)  # Sold between 2000 and 2019
i_ls &= d['p_e_ls'].eq(0)  # Not protected by an easement at the time of sale

# Select rows based on the criteria
X_ls = d.loc[i_ls, XCOLS_LS].copy()
y_ls = d.loc[i_ls, 'ls_price_per_ha_log'].copy()

# Convert 'ls_year' values to integers
X_ls['ls_year'] = X_ls['ls_year'].astype(int)
In [10]:
# Define the selection criteria for X_fc and y_fc
i_fc = d[XCOLS_FC].isnull().sum(1).eq(0)
i_fc &= d['p_prot'].le(5)  # Not protected (p_prot <= 5%)

# Select rows based on the criteria
X_fc = d.loc[i_fc, XCOLS_FC].copy()
y_fc = d.loc[i_fc, 'p_f_change'].copy()
In [11]:
import statsmodels.api as sm
model_ls_ols = sm.OLS(y_ls, sm.add_constant(X_ls)).fit()
model_ls_ols.summary()
Out[11]:
OLS Regression Results
Dep. Variable: ls_price_per_ha_log R-squared: 0.455
Model: OLS Adj. R-squared: 0.455
Method: Least Squares F-statistic: 1687.
Date: Thu, 14 Dec 2023 Prob (F-statistic): 0.00
Time: 17:18:24 Log-Likelihood: -78212.
No. Observations: 44424 AIC: 1.565e+05
Df Residuals: 44401 BIC: 1.567e+05
Df Model: 22
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const -52.9732 2.823 -18.764 0.000 -58.507 -47.440
slope -0.0266 0.002 -10.791 0.000 -0.031 -0.022
travel -0.0006 0.000 -6.079 0.000 -0.001 -0.000
p_wet -0.0050 0.000 -24.097 0.000 -0.005 -0.005
coast_2500 0.0298 0.001 29.846 0.000 0.028 0.032
river_frontage -2.028e-05 7.75e-05 -0.262 0.794 -0.000 0.000
lake_frontage -0.0007 5.09e-05 -12.880 0.000 -0.001 -0.001
pop_dens_tract_2012-2016 -3.011e-05 1.99e-05 -1.511 0.131 -6.92e-05 8.94e-06
hh_inc_avg_tract_2012-2016 1.322e-05 1.88e-07 70.248 0.000 1.29e-05 1.36e-05
p_f_2000 -0.0059 0.000 -16.005 0.000 -0.007 -0.005
p_prot_2000_500 -2.324e-05 0.001 -0.046 0.964 -0.001 0.001
p_d_2000_200 0.0124 0.001 17.975 0.000 0.011 0.014
p_d_2000_1000 -0.0002 0.001 -0.250 0.803 -0.002 0.002
p_d_2000_5000 0.0154 0.001 16.528 0.000 0.014 0.017
p_dev_open 0.0036 0.001 6.977 0.000 0.003 0.005
p_dev_low 0.0075 0.001 9.104 0.000 0.006 0.009
p_dev_medium 0.0150 0.001 20.060 0.000 0.013 0.016
p_dev_high 0.0045 0.001 5.636 0.000 0.003 0.006
p_crops -0.0038 0.001 -4.302 0.000 -0.006 -0.002
p_pasture -0.0002 0.001 -0.382 0.702 -0.001 0.001
p_grassland 0.0024 0.001 2.224 0.026 0.000 0.004
bld_area_per_ha 1.86e-05 6.64e-07 28.024 0.000 1.73e-05 1.99e-05
ls_year 0.0315 0.001 22.385 0.000 0.029 0.034
Omnibus: 7448.307 Durbin-Watson: 1.958
Prob(Omnibus): 0.000 Jarque-Bera (JB): 34364.698
Skew: -0.752 Prob(JB): 0.00
Kurtosis: 7.038 Cond. No. 4.82e+07


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.82e+07. This might indicate that there are
strong multicollinearity or other numerical problems.
In [12]:
fig, ax = plt.subplots(figsize=(17, 10))
dp = shp_c.join(model_ls_ols.fittedvalues.rename('pred'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp.plot(
    'pred',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.axis('off')
Out[12]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [13]:
fig, ax = plt.subplots(figsize=(17, 10))
rp = shp_c.join(model_ls_ols.resid.rename('resid'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
rp.plot(
    'resid',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.axis('off')
Out[13]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [14]:
df = X_ls[['ls_year']].join(model_ls_ols.resid.rename('resid'))
df.groupby('ls_year')['resid'].mean().plot()
Out[14]:
<Axes: xlabel='ls_year'>
In [15]:
from sklearn.metrics import mean_squared_error
mean_squared_error(y_ls, model_ls_ols.fittedvalues)
Out[15]:
1.980329829753318
In [16]:
from sklearn.preprocessing import OneHotEncoder

enc = OneHotEncoder(categories='auto').fit(X_ls[['ls_year']])
dummies_array = enc.transform(X_ls[['ls_year']]).toarray()
dummies_df = pd.DataFrame(
    dummies_array,
    index=X_ls.index,
    columns=[v[8:] for v in enc.get_feature_names_out()],
)
In [17]:
X_ls2 = X_ls.drop(columns='ls_year').join(dummies_df).drop(columns='2000')
In [18]:
model_ls_ols2 = sm.OLS(y_ls, sm.add_constant(X_ls2)).fit()
model_ls_ols2.summary()
Out[18]:
OLS Regression Results
Dep. Variable: ls_price_per_ha_log R-squared: 0.462
Model: OLS Adj. R-squared: 0.461
Method: Least Squares F-statistic: 1002.
Date: Thu, 14 Dec 2023 Prob (F-statistic): 0.00
Time: 17:18:31 Log-Likelihood: -77948.
No. Observations: 44424 AIC: 1.560e+05
Df Residuals: 44385 BIC: 1.563e+05
Df Model: 38
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 9.7614 0.053 182.862 0.000 9.657 9.866
slope -0.0274 0.002 -11.177 0.000 -0.032 -0.023
travel -0.0007 0.000 -6.493 0.000 -0.001 -0.000
p_wet -0.0051 0.000 -24.486 0.000 -0.005 -0.005
coast_2500 0.0300 0.001 30.207 0.000 0.028 0.032
river_frontage -3.626e-05 7.71e-05 -0.471 0.638 -0.000 0.000
lake_frontage -0.0007 5.07e-05 -12.922 0.000 -0.001 -0.001
pop_dens_tract_2012-2016 -3.376e-05 1.98e-05 -1.704 0.088 -7.26e-05 5.08e-06
hh_inc_avg_tract_2012-2016 1.316e-05 1.87e-07 70.301 0.000 1.28e-05 1.35e-05
p_f_2000 -0.0058 0.000 -16.027 0.000 -0.007 -0.005
p_prot_2000_500 8.67e-05 0.001 0.172 0.864 -0.001 0.001
p_d_2000_200 0.0126 0.001 18.356 0.000 0.011 0.014
p_d_2000_1000 -0.0004 0.001 -0.486 0.627 -0.002 0.001
p_d_2000_5000 0.0155 0.001 16.736 0.000 0.014 0.017
p_dev_open 0.0034 0.001 6.633 0.000 0.002 0.004
p_dev_low 0.0073 0.001 9.010 0.000 0.006 0.009
p_dev_medium 0.0146 0.001 19.662 0.000 0.013 0.016
p_dev_high 0.0040 0.001 5.138 0.000 0.002 0.006
p_crops -0.0039 0.001 -4.439 0.000 -0.006 -0.002
p_pasture -0.0001 0.001 -0.253 0.801 -0.001 0.001
p_grassland 0.0021 0.001 2.037 0.042 8.14e-05 0.004
bld_area_per_ha 1.866e-05 6.6e-07 28.272 0.000 1.74e-05 2e-05
2001 0.0837 0.039 2.160 0.031 0.008 0.160
2002 0.0817 0.038 2.152 0.031 0.007 0.156
2003 0.2283 0.038 6.033 0.000 0.154 0.302
2004 0.4424 0.037 12.035 0.000 0.370 0.514
2005 0.6535 0.037 17.885 0.000 0.582 0.725
2006 0.6194 0.037 16.565 0.000 0.546 0.693
2007 0.6743 0.038 17.898 0.000 0.600 0.748
2008 0.5678 0.039 14.575 0.000 0.491 0.644
2009 0.5697 0.040 14.301 0.000 0.492 0.648
2010 0.5044 0.039 13.089 0.000 0.429 0.580
2011 0.5055 0.040 12.682 0.000 0.427 0.584
2012 0.5223 0.039 13.351 0.000 0.446 0.599
2013 0.5588 0.039 14.379 0.000 0.483 0.635
2014 0.5453 0.039 13.931 0.000 0.469 0.622
2015 0.6266 0.041 15.312 0.000 0.546 0.707
2016 0.6110 0.043 14.099 0.000 0.526 0.696
2017 0.7002 0.109 6.409 0.000 0.486 0.914
Omnibus: 7598.673 Durbin-Watson: 1.954
Prob(Omnibus): 0.000 Jarque-Bera (JB): 37366.196
Skew: -0.749 Prob(JB): 0.00
Kurtosis: 7.236 Cond. No. 2.16e+06


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.16e+06. This might indicate that there are
strong multicollinearity or other numerical problems.
In [19]:
mean_squared_error(y_ls, model_ls_ols2.fittedvalues)
Out[19]:
1.9570070916504736
In [20]:
model_ls_ols2.params.loc['2001':'2017'].plot()
Out[20]:
<Axes: >
In [21]:
from itertools import combinations
X_ls3 = X_ls.drop(columns='ls_year').copy()
cols = X_ls3.columns
for v in cols:
    X_ls3[v + '^2'] = X_ls3[v] ** 2
for v1, v2 in combinations(cols,2):
    X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
X_ls3 = X_ls3.join(dummies_df)
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
In [22]:
model_ls_ols3 = sm.OLS(y_ls, sm.add_constant(X_ls3)).fit()
model_ls_ols3.summary()
Out[22]:
OLS Regression Results
Dep. Variable: ls_price_per_ha_log R-squared: 0.537
Model: OLS Adj. R-squared: 0.534
Method: Least Squares F-statistic: 190.2
Date: Thu, 14 Dec 2023 Prob (F-statistic): 0.00
Time: 17:18:33 Log-Likelihood: -74613.
No. Observations: 44424 AIC: 1.498e+05
Df Residuals: 44154 BIC: 1.521e+05
Df Model: 269
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 9.6012 0.191 50.254 0.000 9.227 9.976
slope -0.0174 0.020 -0.877 0.380 -0.056 0.021
travel -0.0016 0.001 -1.762 0.078 -0.003 0.000
p_wet -0.0138 0.001 -10.100 0.000 -0.016 -0.011
coast_2500 0.0679 0.009 7.977 0.000 0.051 0.085
river_frontage 0.0053 0.001 9.006 0.000 0.004 0.006
lake_frontage -0.0021 0.000 -4.445 0.000 -0.003 -0.001
pop_dens_tract_2012-2016 0.0009 0.000 4.706 0.000 0.001 0.001
hh_inc_avg_tract_2012-2016 2.146e-05 1.42e-06 15.101 0.000 1.87e-05 2.42e-05
p_f_2000 -0.0146 0.003 -5.670 0.000 -0.020 -0.010
p_prot_2000_500 0.0017 0.004 0.484 0.629 -0.005 0.009
p_d_2000_200 0.0138 0.004 3.101 0.002 0.005 0.022
p_d_2000_1000 0.0258 0.006 4.535 0.000 0.015 0.037
p_d_2000_5000 -0.0066 0.006 -1.053 0.292 -0.019 0.006
p_dev_open 0.0025 0.003 0.817 0.414 -0.003 0.008
p_dev_low -0.0023 0.005 -0.430 0.667 -0.013 0.008
p_dev_medium 0.0095 0.006 1.580 0.114 -0.002 0.021
p_dev_high -0.0162 0.008 -1.902 0.057 -0.033 0.000
p_crops 0.0026 0.007 0.348 0.728 -0.012 0.017
p_pasture -0.0241 0.004 -6.016 0.000 -0.032 -0.016
p_grassland -0.0202 0.007 -2.732 0.006 -0.035 -0.006
bld_area_per_ha 5.731e-05 5.32e-06 10.776 0.000 4.69e-05 6.77e-05
slope^2 0.0009 0.000 2.055 0.040 4.25e-05 0.002
travel^2 -2.908e-06 9.37e-07 -3.104 0.002 -4.74e-06 -1.07e-06
p_wet^2 1.542e-05 4.45e-06 3.466 0.001 6.7e-06 2.41e-05
coast_2500^2 -0.0002 7.04e-05 -3.222 0.001 -0.000 -8.88e-05
river_frontage^2 -2.876e-07 5.29e-08 -5.439 0.000 -3.91e-07 -1.84e-07
lake_frontage^2 7.263e-08 2.17e-08 3.354 0.001 3.02e-08 1.15e-07
pop_dens_tract_2012-2016^2 -1.153e-08 1.01e-08 -1.140 0.254 -3.14e-08 8.29e-09
hh_inc_avg_tract_2012-2016^2 -2.954e-11 2.85e-12 -10.383 0.000 -3.51e-11 -2.4e-11
p_f_2000^2 6.417e-06 1.33e-05 0.481 0.631 -1.97e-05 3.26e-05
p_prot_2000_500^2 -2.069e-05 2.16e-05 -0.956 0.339 -6.31e-05 2.17e-05
p_d_2000_200^2 -1.46e-05 4.05e-05 -0.361 0.718 -9.39e-05 6.47e-05
p_d_2000_1000^2 4.381e-05 6.27e-05 0.699 0.485 -7.9e-05 0.000
p_d_2000_5000^2 -9.542e-05 7.49e-05 -1.274 0.203 -0.000 5.14e-05
p_dev_open^2 -3.831e-05 1.97e-05 -1.949 0.051 -7.68e-05 2.22e-07
p_dev_low^2 -1.355e-05 4.15e-05 -0.326 0.744 -9.5e-05 6.79e-05
p_dev_medium^2 -0.0001 4.96e-05 -2.144 0.032 -0.000 -9.13e-06
p_dev_high^2 0.0002 7.41e-05 2.203 0.028 1.8e-05 0.000
p_crops^2 4.586e-05 5.49e-05 0.836 0.403 -6.17e-05 0.000
p_pasture^2 0.0002 2.78e-05 7.047 0.000 0.000 0.000
p_grassland^2 0.0001 4.84e-05 2.277 0.023 1.54e-05 0.000
bld_area_per_ha^2 -8.609e-11 6.27e-12 -13.733 0.000 -9.84e-11 -7.38e-11
slope*travel 8.351e-05 3.11e-05 2.688 0.007 2.26e-05 0.000
slope*p_wet 0.0006 0.000 4.816 0.000 0.000 0.001
slope*coast_2500 0.0010 0.001 1.923 0.054 -1.93e-05 0.002
slope*river_frontage 0.0001 2.96e-05 3.935 0.000 5.84e-05 0.000
slope*lake_frontage 6.781e-05 3.13e-05 2.164 0.030 6.39e-06 0.000
slope*pop_dens_tract_2012-2016 -1.725e-05 9.02e-06 -1.911 0.056 -3.49e-05 4.4e-07
slope*hh_inc_avg_tract_2012-2016 2.845e-07 8.05e-08 3.536 0.000 1.27e-07 4.42e-07
slope*p_f_2000 -0.0008 0.000 -4.936 0.000 -0.001 -0.000
slope*p_prot_2000_500 -0.0002 0.000 -1.440 0.150 -0.001 8.32e-05
slope*p_d_2000_200 0.0005 0.000 1.549 0.121 -0.000 0.001
slope*p_d_2000_1000 0.0001 0.000 0.365 0.715 -0.001 0.001
slope*p_d_2000_5000 0.0016 0.000 4.305 0.000 0.001 0.002
slope*p_dev_open -0.0003 0.000 -1.820 0.069 -0.001 2.67e-05
slope*p_dev_low -0.0006 0.000 -1.556 0.120 -0.001 0.000
slope*p_dev_medium -0.0011 0.000 -2.975 0.003 -0.002 -0.000
slope*p_dev_high -0.0018 0.000 -4.194 0.000 -0.003 -0.001
slope*p_crops -0.0001 0.000 -0.275 0.783 -0.001 0.001
slope*p_pasture -0.0003 0.000 -1.343 0.179 -0.001 0.000
slope*p_grassland 0.0003 0.001 0.488 0.626 -0.001 0.001
slope*bld_area_per_ha 1.225e-06 3.13e-07 3.911 0.000 6.11e-07 1.84e-06
travel*p_wet -1.254e-07 4.01e-06 -0.031 0.975 -7.98e-06 7.73e-06
travel*coast_2500 3.707e-05 1.91e-05 1.941 0.052 -3.65e-07 7.45e-05
travel*river_frontage -2.083e-06 2.43e-06 -0.857 0.391 -6.85e-06 2.68e-06
travel*lake_frontage 3.203e-06 1.15e-06 2.793 0.005 9.56e-07 5.45e-06
travel*pop_dens_tract_2012-2016 -1.191e-05 1.46e-06 -8.140 0.000 -1.48e-05 -9.04e-06
travel*hh_inc_avg_tract_2012-2016 -6.461e-09 5.47e-09 -1.182 0.237 -1.72e-08 4.25e-09
travel*p_f_2000 1.716e-05 6.69e-06 2.563 0.010 4.04e-06 3.03e-05
travel*p_prot_2000_500 -7.408e-06 6.99e-06 -1.060 0.289 -2.11e-05 6.29e-06
travel*p_d_2000_200 3.81e-05 1.42e-05 2.685 0.007 1.03e-05 6.59e-05
travel*p_d_2000_1000 -2.652e-05 2.22e-05 -1.197 0.231 -6.99e-05 1.69e-05
travel*p_d_2000_5000 0.0002 2.94e-05 7.924 0.000 0.000 0.000
travel*p_dev_open 4.246e-05 9.85e-06 4.309 0.000 2.31e-05 6.18e-05
travel*p_dev_low -1.953e-05 2.28e-05 -0.856 0.392 -6.42e-05 2.52e-05
travel*p_dev_medium -0.0002 3.81e-05 -5.130 0.000 -0.000 -0.000
travel*p_dev_high -0.0002 5.85e-05 -3.387 0.001 -0.000 -8.35e-05
travel*p_crops -4.288e-05 2.35e-05 -1.821 0.069 -8.9e-05 3.27e-06
travel*p_pasture 2.318e-05 8.07e-06 2.871 0.004 7.36e-06 3.9e-05
travel*p_grassland 4.614e-05 1.66e-05 2.778 0.005 1.36e-05 7.87e-05
travel*bld_area_per_ha 6.117e-07 4.17e-08 14.682 0.000 5.3e-07 6.93e-07
p_wet*coast_2500 -0.0001 2.58e-05 -4.543 0.000 -0.000 -6.65e-05
p_wet*river_frontage -8.636e-06 2.4e-06 -3.598 0.000 -1.33e-05 -3.93e-06
p_wet*lake_frontage 4.736e-06 1.98e-06 2.397 0.017 8.64e-07 8.61e-06
p_wet*pop_dens_tract_2012-2016 1.93e-06 8.17e-07 2.364 0.018 3.3e-07 3.53e-06
p_wet*hh_inc_avg_tract_2012-2016 1.195e-08 5.26e-09 2.272 0.023 1.64e-09 2.23e-08
p_wet*p_f_2000 3.273e-05 8.24e-06 3.974 0.000 1.66e-05 4.89e-05
p_wet*p_prot_2000_500 -2.008e-05 1.45e-05 -1.381 0.167 -4.86e-05 8.42e-06
p_wet*p_d_2000_200 0.0001 1.9e-05 6.371 0.000 8.38e-05 0.000
p_wet*p_d_2000_1000 -2.509e-05 2.44e-05 -1.028 0.304 -7.29e-05 2.28e-05
p_wet*p_d_2000_5000 -4.172e-06 2.64e-05 -0.158 0.874 -5.59e-05 4.76e-05
p_wet*p_dev_open -1.866e-05 1.36e-05 -1.369 0.171 -4.54e-05 8.06e-06
p_wet*p_dev_low 5.079e-05 2.48e-05 2.051 0.040 2.26e-06 9.93e-05
p_wet*p_dev_medium 2.107e-05 3.17e-05 0.665 0.506 -4.1e-05 8.31e-05
p_wet*p_dev_high -4.848e-05 5.09e-05 -0.953 0.341 -0.000 5.12e-05
p_wet*p_crops 6.887e-05 2.21e-05 3.121 0.002 2.56e-05 0.000
p_wet*p_pasture 9.225e-05 2.37e-05 3.891 0.000 4.58e-05 0.000
p_wet*p_grassland 7.173e-05 4.21e-05 1.704 0.088 -1.08e-05 0.000
p_wet*bld_area_per_ha 2.256e-07 4.47e-08 5.051 0.000 1.38e-07 3.13e-07
coast_2500*river_frontage -2.147e-05 1.13e-05 -1.895 0.058 -4.37e-05 7.37e-07
coast_2500*lake_frontage 1.275e-05 7.43e-06 1.717 0.086 -1.8e-06 2.73e-05
coast_2500*pop_dens_tract_2012-2016 -3.196e-06 2.7e-06 -1.184 0.236 -8.49e-06 2.09e-06
coast_2500*hh_inc_avg_tract_2012-2016 -2.148e-07 3.71e-08 -5.783 0.000 -2.88e-07 -1.42e-07
coast_2500*p_f_2000 -1.858e-05 3.99e-05 -0.465 0.642 -9.69e-05 5.97e-05
coast_2500*p_prot_2000_500 -6.632e-05 5.99e-05 -1.107 0.268 -0.000 5.11e-05
coast_2500*p_d_2000_200 -5.542e-05 7.64e-05 -0.725 0.468 -0.000 9.43e-05
coast_2500*p_d_2000_1000 -0.0002 0.000 -1.516 0.129 -0.000 4.61e-05
coast_2500*p_d_2000_5000 -0.0002 0.000 -1.086 0.277 -0.000 0.000
coast_2500*p_dev_open -0.0001 6.39e-05 -1.761 0.078 -0.000 1.27e-05
coast_2500*p_dev_low -7.24e-05 9.53e-05 -0.760 0.448 -0.000 0.000
coast_2500*p_dev_medium 6.551e-05 0.000 0.653 0.514 -0.000 0.000
coast_2500*p_dev_high -8.371e-05 0.000 -0.800 0.424 -0.000 0.000
coast_2500*p_crops 0.0001 0.000 0.495 0.621 -0.000 0.001
coast_2500*p_pasture -0.0003 0.000 -2.069 0.039 -0.001 -1.54e-05
coast_2500*p_grassland 8.18e-05 7.2e-05 1.137 0.256 -5.93e-05 0.000
coast_2500*bld_area_per_ha -2.167e-07 8.71e-08 -2.489 0.013 -3.87e-07 -4.61e-08
river_frontage*lake_frontage 2.478e-07 2.08e-07 1.189 0.235 -1.61e-07 6.56e-07
river_frontage*pop_dens_tract_2012-2016 -1.523e-07 2.72e-07 -0.560 0.576 -6.85e-07 3.81e-07
river_frontage*hh_inc_avg_tract_2012-2016 -1.595e-08 3.13e-09 -5.090 0.000 -2.21e-08 -9.81e-09
river_frontage*p_f_2000 -3.541e-05 4.67e-06 -7.575 0.000 -4.46e-05 -2.62e-05
river_frontage*p_prot_2000_500 -3.768e-06 7.39e-06 -0.510 0.610 -1.83e-05 1.07e-05
river_frontage*p_d_2000_200 -4.3e-05 1.09e-05 -3.953 0.000 -6.43e-05 -2.17e-05
river_frontage*p_d_2000_1000 -1.34e-05 1.06e-05 -1.268 0.205 -3.41e-05 7.31e-06
river_frontage*p_d_2000_5000 1.521e-05 1.08e-05 1.411 0.158 -5.92e-06 3.63e-05
river_frontage*p_dev_open -1.759e-05 5.73e-06 -3.072 0.002 -2.88e-05 -6.37e-06
river_frontage*p_dev_low -7.647e-06 1.68e-05 -0.456 0.648 -4.05e-05 2.52e-05
river_frontage*p_dev_medium 1.223e-05 1.82e-05 0.672 0.502 -2.35e-05 4.79e-05
river_frontage*p_dev_high -3.513e-05 1.5e-05 -2.339 0.019 -6.46e-05 -5.7e-06
river_frontage*p_crops 1.597e-05 8.8e-06 1.815 0.070 -1.28e-06 3.32e-05
river_frontage*p_pasture -5.539e-05 8.58e-06 -6.453 0.000 -7.22e-05 -3.86e-05
river_frontage*p_grassland -5.353e-05 1.64e-05 -3.271 0.001 -8.56e-05 -2.15e-05
river_frontage*bld_area_per_ha -4.176e-08 1.28e-08 -3.269 0.001 -6.68e-08 -1.67e-08
lake_frontage*pop_dens_tract_2012-2016 -3.489e-07 2.48e-07 -1.407 0.160 -8.35e-07 1.37e-07
lake_frontage*hh_inc_avg_tract_2012-2016 -2.582e-10 1.73e-09 -0.149 0.881 -3.65e-09 3.13e-09
lake_frontage*p_f_2000 -1.533e-06 3.83e-06 -0.400 0.689 -9.03e-06 5.97e-06
lake_frontage*p_prot_2000_500 1.11e-06 3.58e-06 0.310 0.757 -5.91e-06 8.13e-06
lake_frontage*p_d_2000_200 1.067e-05 7.66e-06 1.393 0.164 -4.35e-06 2.57e-05
lake_frontage*p_d_2000_1000 -7.129e-06 8.98e-06 -0.794 0.427 -2.47e-05 1.05e-05
lake_frontage*p_d_2000_5000 3.146e-05 8.46e-06 3.719 0.000 1.49e-05 4.8e-05
lake_frontage*p_dev_open -2.499e-07 4.82e-06 -0.052 0.959 -9.71e-06 9.21e-06
lake_frontage*p_dev_low 1.356e-05 1.04e-05 1.302 0.193 -6.85e-06 3.4e-05
lake_frontage*p_dev_medium 1.452e-05 1.12e-05 1.295 0.195 -7.45e-06 3.65e-05
lake_frontage*p_dev_high 6.245e-06 9.65e-06 0.647 0.517 -1.27e-05 2.52e-05
lake_frontage*p_crops -1.792e-05 7.24e-06 -2.476 0.013 -3.21e-05 -3.73e-06
lake_frontage*p_pasture -3.026e-06 6.68e-06 -0.453 0.651 -1.61e-05 1.01e-05
lake_frontage*p_grassland -6.274e-06 1.12e-05 -0.562 0.574 -2.81e-05 1.56e-05
lake_frontage*bld_area_per_ha -2.614e-08 1.04e-08 -2.523 0.012 -4.65e-08 -5.83e-09
pop_dens_tract_2012-2016*hh_inc_avg_tract_2012-2016 2.901e-09 7.01e-10 4.140 0.000 1.53e-09 4.27e-09
pop_dens_tract_2012-2016*p_f_2000 -1.748e-06 1.56e-06 -1.120 0.263 -4.81e-06 1.31e-06
pop_dens_tract_2012-2016*p_prot_2000_500 3.828e-08 2.34e-06 0.016 0.987 -4.54e-06 4.62e-06
pop_dens_tract_2012-2016*p_d_2000_200 3.151e-06 1.87e-06 1.686 0.092 -5.13e-07 6.81e-06
pop_dens_tract_2012-2016*p_d_2000_1000 -8.531e-06 2.02e-06 -4.232 0.000 -1.25e-05 -4.58e-06
pop_dens_tract_2012-2016*p_d_2000_5000 -2.691e-06 2.03e-06 -1.326 0.185 -6.67e-06 1.29e-06
pop_dens_tract_2012-2016*p_dev_open 1.069e-06 1.89e-06 0.565 0.572 -2.64e-06 4.78e-06
pop_dens_tract_2012-2016*p_dev_low -5.17e-06 2.22e-06 -2.328 0.020 -9.52e-06 -8.18e-07
pop_dens_tract_2012-2016*p_dev_medium -1.27e-06 1.67e-06 -0.760 0.447 -4.55e-06 2.01e-06
pop_dens_tract_2012-2016*p_dev_high -3.519e-06 1.61e-06 -2.181 0.029 -6.68e-06 -3.57e-07
pop_dens_tract_2012-2016*p_crops -7.678e-06 8.48e-06 -0.905 0.365 -2.43e-05 8.95e-06
pop_dens_tract_2012-2016*p_pasture 6.198e-06 4.12e-06 1.506 0.132 -1.87e-06 1.43e-05
pop_dens_tract_2012-2016*p_grassland 3.62e-06 5.36e-06 0.676 0.499 -6.88e-06 1.41e-05
pop_dens_tract_2012-2016*bld_area_per_ha 1.131e-09 4.81e-10 2.354 0.019 1.89e-10 2.07e-09
hh_inc_avg_tract_2012-2016*p_f_2000 1.149e-08 9.55e-09 1.204 0.229 -7.22e-09 3.02e-08
hh_inc_avg_tract_2012-2016*p_prot_2000_500 -2.862e-08 1.29e-08 -2.213 0.027 -5.4e-08 -3.27e-09
hh_inc_avg_tract_2012-2016*p_d_2000_200 -9.507e-08 1.84e-08 -5.168 0.000 -1.31e-07 -5.9e-08
hh_inc_avg_tract_2012-2016*p_d_2000_1000 -4.277e-08 2.38e-08 -1.796 0.073 -8.94e-08 3.91e-09
hh_inc_avg_tract_2012-2016*p_d_2000_5000 -4.514e-08 2.39e-08 -1.889 0.059 -9.2e-08 1.7e-09
hh_inc_avg_tract_2012-2016*p_dev_open -1.055e-08 1.25e-08 -0.847 0.397 -3.5e-08 1.39e-08
hh_inc_avg_tract_2012-2016*p_dev_low -7.634e-09 1.98e-08 -0.386 0.699 -4.64e-08 3.11e-08
hh_inc_avg_tract_2012-2016*p_dev_medium -3.384e-08 1.98e-08 -1.706 0.088 -7.27e-08 5.04e-09
hh_inc_avg_tract_2012-2016*p_dev_high 1.765e-09 2.32e-08 0.076 0.939 -4.37e-08 4.72e-08
hh_inc_avg_tract_2012-2016*p_crops -1.002e-07 4.04e-08 -2.480 0.013 -1.79e-07 -2.1e-08
hh_inc_avg_tract_2012-2016*p_pasture -7.725e-09 1.49e-08 -0.519 0.604 -3.69e-08 2.14e-08
hh_inc_avg_tract_2012-2016*p_grassland -2.71e-10 3.67e-08 -0.007 0.994 -7.21e-08 7.16e-08
hh_inc_avg_tract_2012-2016*bld_area_per_ha 7.023e-11 1.58e-11 4.448 0.000 3.93e-11 1.01e-10
p_f_2000*p_prot_2000_500 3.211e-05 2.64e-05 1.218 0.223 -1.96e-05 8.38e-05
p_f_2000*p_d_2000_200 9.3e-05 3.2e-05 2.906 0.004 3.03e-05 0.000
p_f_2000*p_d_2000_1000 -0.0001 4.13e-05 -3.200 0.001 -0.000 -5.12e-05
p_f_2000*p_d_2000_5000 0.0001 4.61e-05 3.159 0.002 5.53e-05 0.000
p_f_2000*p_dev_open 0.0001 2.03e-05 5.591 0.000 7.38e-05 0.000
p_f_2000*p_dev_low 0.0002 3.46e-05 6.589 0.000 0.000 0.000
p_f_2000*p_dev_medium -0.0001 4.31e-05 -2.664 0.008 -0.000 -3.04e-05
p_f_2000*p_dev_high -3.146e-05 5.37e-05 -0.586 0.558 -0.000 7.38e-05
p_f_2000*p_crops 0.0002 4.98e-05 3.062 0.002 5.49e-05 0.000
p_f_2000*p_pasture 0.0002 3e-05 5.381 0.000 0.000 0.000
p_f_2000*p_grassland 8.42e-05 4.11e-05 2.050 0.040 3.69e-06 0.000
p_f_2000*bld_area_per_ha 7.959e-07 5.17e-08 15.393 0.000 6.95e-07 8.97e-07
p_prot_2000_500*p_d_2000_200 0.0001 5.55e-05 2.055 0.040 5.28e-06 0.000
p_prot_2000_500*p_d_2000_1000 -0.0002 6.83e-05 -3.428 0.001 -0.000 -0.000
p_prot_2000_500*p_d_2000_5000 0.0001 6.94e-05 2.137 0.033 1.23e-05 0.000
p_prot_2000_500*p_dev_open 3.756e-05 3.99e-05 0.942 0.346 -4.06e-05 0.000
p_prot_2000_500*p_dev_low -6.002e-05 6.37e-05 -0.942 0.346 -0.000 6.48e-05
p_prot_2000_500*p_dev_medium 0.0002 6.75e-05 3.600 0.000 0.000 0.000
p_prot_2000_500*p_dev_high 0.0001 8.26e-05 1.396 0.163 -4.66e-05 0.000
p_prot_2000_500*p_crops -8.785e-05 6.83e-05 -1.286 0.198 -0.000 4.6e-05
p_prot_2000_500*p_pasture 0.0001 3.71e-05 3.228 0.001 4.71e-05 0.000
p_prot_2000_500*p_grassland -1.474e-06 7.19e-05 -0.021 0.984 -0.000 0.000
p_prot_2000_500*bld_area_per_ha -2.722e-07 6.34e-08 -4.293 0.000 -3.97e-07 -1.48e-07
p_d_2000_200*p_d_2000_1000 -0.0003 7.95e-05 -3.950 0.000 -0.000 -0.000
p_d_2000_200*p_d_2000_5000 9.331e-05 7.77e-05 1.201 0.230 -5.9e-05 0.000
p_d_2000_200*p_dev_open -0.0001 4.07e-05 -3.487 0.000 -0.000 -6.22e-05
p_d_2000_200*p_dev_low 0.0001 5.94e-05 1.762 0.078 -1.17e-05 0.000
p_d_2000_200*p_dev_medium 0.0002 5.95e-05 3.596 0.000 9.73e-05 0.000
p_d_2000_200*p_dev_high 0.0002 6.4e-05 2.538 0.011 3.7e-05 0.000
p_d_2000_200*p_crops -9.669e-05 8.61e-05 -1.123 0.261 -0.000 7.21e-05
p_d_2000_200*p_pasture 0.0001 4.83e-05 2.984 0.003 4.95e-05 0.000
p_d_2000_200*p_grassland 3.766e-05 7.54e-05 0.500 0.617 -0.000 0.000
p_d_2000_200*bld_area_per_ha -1.715e-07 5.94e-08 -2.888 0.004 -2.88e-07 -5.51e-08
p_d_2000_1000*p_d_2000_5000 -3.977e-05 0.000 -0.381 0.703 -0.000 0.000
p_d_2000_1000*p_dev_open -3.03e-05 5.62e-05 -0.539 0.590 -0.000 7.99e-05
p_d_2000_1000*p_dev_low -1.616e-05 8.09e-05 -0.200 0.842 -0.000 0.000
p_d_2000_1000*p_dev_medium 0.0002 7.04e-05 3.344 0.001 9.74e-05 0.000
p_d_2000_1000*p_dev_high 0.0001 7.28e-05 1.439 0.150 -3.79e-05 0.000
p_d_2000_1000*p_crops 0.0001 0.000 1.055 0.292 -0.000 0.000
p_d_2000_1000*p_pasture -0.0002 7.03e-05 -2.595 0.009 -0.000 -4.46e-05
p_d_2000_1000*p_grassland -3.903e-05 0.000 -0.345 0.730 -0.000 0.000
p_d_2000_1000*bld_area_per_ha -2.973e-07 6.26e-08 -4.746 0.000 -4.2e-07 -1.75e-07
p_d_2000_5000*p_dev_open 5.358e-05 6.31e-05 0.849 0.396 -7.01e-05 0.000
p_d_2000_5000*p_dev_low 3.066e-05 9.16e-05 0.335 0.738 -0.000 0.000
p_d_2000_5000*p_dev_medium -3.632e-05 7.72e-05 -0.470 0.638 -0.000 0.000
p_d_2000_5000*p_dev_high 9.372e-05 7.82e-05 1.199 0.231 -5.95e-05 0.000
p_d_2000_5000*p_crops -0.0003 0.000 -2.153 0.031 -0.001 -2.63e-05
p_d_2000_5000*p_pasture -3.688e-05 8.66e-05 -0.426 0.670 -0.000 0.000
p_d_2000_5000*p_grassland 0.0002 0.000 1.210 0.226 -0.000 0.000
p_d_2000_5000*bld_area_per_ha 4.189e-07 6.01e-08 6.969 0.000 3.01e-07 5.37e-07
p_dev_open*p_dev_low 1.793e-05 5.24e-05 0.342 0.732 -8.47e-05 0.000
p_dev_open*p_dev_medium 0.0001 7.96e-05 1.872 0.061 -7.02e-06 0.000
p_dev_open*p_dev_high 0.0002 0.000 1.689 0.091 -3.7e-05 0.000
p_dev_open*p_crops -0.0002 0.000 -1.399 0.162 -0.001 8.53e-05
p_dev_open*p_pasture -2.518e-05 5.84e-05 -0.431 0.666 -0.000 8.92e-05
p_dev_open*p_grassland -5.465e-06 0.000 -0.049 0.961 -0.000 0.000
p_dev_open*bld_area_per_ha -4.416e-07 9.98e-08 -4.425 0.000 -6.37e-07 -2.46e-07
p_dev_low*p_dev_medium 0.0001 6.87e-05 1.768 0.077 -1.32e-05 0.000
p_dev_low*p_dev_high 4.714e-05 0.000 0.395 0.693 -0.000 0.000
p_dev_low*p_crops 0.0002 0.000 1.056 0.291 -0.000 0.001
p_dev_low*p_pasture 0.0002 8.25e-05 2.157 0.031 1.63e-05 0.000
p_dev_low*p_grassland 0.0001 0.000 0.787 0.431 -0.000 0.000
p_dev_low*bld_area_per_ha -5.213e-07 6.86e-08 -7.604 0.000 -6.56e-07 -3.87e-07
p_dev_medium*p_dev_high -4.418e-05 7.77e-05 -0.569 0.569 -0.000 0.000
p_dev_medium*p_crops -0.0003 0.000 -1.449 0.147 -0.001 0.000
p_dev_medium*p_pasture -3.13e-05 0.000 -0.219 0.827 -0.000 0.000
p_dev_medium*p_grassland 1.027e-05 0.000 0.049 0.961 -0.000 0.000
p_dev_medium*bld_area_per_ha -4.588e-07 4.47e-08 -10.269 0.000 -5.46e-07 -3.71e-07
p_dev_high*p_crops 9.473e-05 0.000 0.213 0.831 -0.001 0.001
p_dev_high*p_pasture 0.0004 0.000 1.366 0.172 -0.000 0.001
p_dev_high*p_grassland 0.0008 0.000 2.157 0.031 7.59e-05 0.002
p_dev_high*bld_area_per_ha -2.93e-07 4.28e-08 -6.840 0.000 -3.77e-07 -2.09e-07
p_crops*p_pasture 0.0002 0.000 1.755 0.079 -2.3e-05 0.000
p_crops*p_grassland 0.0002 0.000 0.706 0.480 -0.000 0.001
p_crops*bld_area_per_ha 8.229e-07 2.53e-07 3.248 0.001 3.26e-07 1.32e-06
p_pasture*p_grassland 0.0002 0.000 1.445 0.148 -6.36e-05 0.000
p_pasture*bld_area_per_ha 2.303e-07 1.36e-07 1.691 0.091 -3.67e-08 4.97e-07
p_grassland*bld_area_per_ha -6.9e-07 2.36e-07 -2.927 0.003 -1.15e-06 -2.28e-07
2000 0.0681 0.027 2.495 0.013 0.015 0.122
2001 0.1649 0.028 5.944 0.000 0.111 0.219
2002 0.1700 0.027 6.353 0.000 0.118 0.222
2003 0.3136 0.027 11.768 0.000 0.261 0.366
2004 0.5303 0.025 20.807 0.000 0.480 0.580
2005 0.7295 0.025 28.882 0.000 0.680 0.779
2006 0.6817 0.026 26.060 0.000 0.630 0.733
2007 0.7235 0.027 27.221 0.000 0.671 0.776
2008 0.6497 0.028 23.296 0.000 0.595 0.704
2009 0.6333 0.029 21.999 0.000 0.577 0.690
2010 0.5561 0.028 20.218 0.000 0.502 0.610
2011 0.5657 0.029 19.493 0.000 0.509 0.623
2012 0.5587 0.028 19.826 0.000 0.503 0.614
2013 0.5877 0.028 21.097 0.000 0.533 0.642
2014 0.6032 0.028 21.488 0.000 0.548 0.658
2015 0.6792 0.030 22.646 0.000 0.620 0.738
2016 0.6543 0.032 20.180 0.000 0.591 0.718
2017 0.7318 0.094 7.769 0.000 0.547 0.916
Omnibus: 8803.431 Durbin-Watson: 1.979
Prob(Omnibus): 0.000 Jarque-Bera (JB): 59018.463
Skew: -0.788 Prob(JB): 0.00
Kurtosis: 8.422 Cond. No. 2.31e+15


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.31e+15. This might indicate that there are
strong multicollinearity or other numerical problems.
In [23]:
mean_squared_error(y_ls, model_ls_ols3.fittedvalues)
Out[23]:
1.6841173295590273
In [24]:
fig, ax = plt.subplots(figsize=(17, 10))
rp3 = shp_c.join(model_ls_ols3.resid.rename('resid'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
rp3.plot(
    'resid',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Residuals of Land Sale Prices Estimation in Massachusetts", fontsize=18)

ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_ols3_resid.png", dpi=200)
Out[24]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [25]:
fig, ax = plt.subplots(figsize=(17, 10))
dp3 = shp_c.join(model_ls_ols3.fittedvalues.rename('pred'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp3.plot(
    'pred',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Values of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_ols3_pred.png", dpi=200)
Out[25]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [26]:
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import (
    ExtraTreesRegressor, 
    GradientBoostingRegressor, 
    RandomForestRegressor
)
from sklearn.model_selection import train_test_split, KFold
In [27]:
model = LinearRegression()
model.fit(X_ls, y_ls)
Out[27]:
LinearRegression()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
LinearRegression()
In [28]:
model.predict(X_ls)
mean_squared_error(y_ls, model.predict(X_ls))
Out[28]:
1.980329829753318
In [29]:
X_train, X_test, y_train, y_test = train_test_split(X_ls, y_ls, test_size=0.5)
model = LinearRegression()
model.fit(X_train, y_train)
mean_squared_error(y_test, model.predict(X_test))
Out[29]:
2.011834891126498
In [30]:
X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = LinearRegression()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 1.9799
Avg MSE (out-of-sample): 1.984
In [31]:
X = X_ls2
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = LinearRegression()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 1.9565
Avg MSE (out-of-sample): 1.9621
In [32]:
X = X_ls3
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = LinearRegression()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 1.6781
Avg MSE (out-of-sample): 1.7686
In [33]:
X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = GradientBoostingRegressor()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 1.4143
Avg MSE (out-of-sample): 1.4732
In [34]:
X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = RandomForestRegressor(n_estimators=20)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 0.2256
Avg MSE (out-of-sample): 1.4086
In [35]:
X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(n_estimators=50)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 0.0
Avg MSE (out-of-sample): 1.3451
In [36]:
feat_imp = pd.Series(model.feature_importances_, index=X.columns)
feat_imp.sort_values().plot(kind='barh')
Out[36]:
<Axes: >
In [37]:
fig, ax = plt.subplots(figsize=(17, 10))
y_pred = pd.Series(model.predict(X), index=X.index)
resid = (y - y_pred).rename('resid')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p = shp_c.join(resid, how='inner')
shp_p.plot(
    'resid',
    ax=ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)
Out[37]:
<Axes: >
In [38]:
model = ExtraTreesRegressor(n_estimators=50, bootstrap=True, oob_score=True)
model.fit(X, y)
y_pred = pd.Series(model.oob_prediction_, index=X.index)
resid = (y - y_pred).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)
Out[38]:
<Axes: >
In [39]:
model = RandomForestRegressor(n_estimators=50, bootstrap=True, oob_score=True)
model.fit(X, y)
y_pred = pd.Series(model.oob_prediction_, index=X.index)
resid = (y - y_pred).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)
Out[39]:
<Axes: >
In [40]:
shp_c['x'] = shp_c['geometry'].x
shp_c['y'] = shp_c['geometry'].y
coords = shp_c[['x', 'y']]
In [41]:
X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = GradientBoostingRegressor()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 1.37
Avg MSE (out-of-sample): 1.4342
In [42]:
X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = RandomForestRegressor(n_estimators=20)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 0.2167
Avg MSE (out-of-sample): 1.3695
In [43]:
X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(n_estimators=50)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 0.0
Avg MSE (out-of-sample): 1.283
In [44]:
X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = GradientBoostingRegressor(n_estimators=1000)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 1.0249
Avg MSE (out-of-sample): 1.3143
In [45]:
X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = RandomForestRegressor(n_estimators=100, oob_score=True)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 0.1833
Avg MSE (out-of-sample): 1.3001
In [46]:
X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(
        n_estimators=100, bootstrap=True, oob_score=True
    )
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 0.1819
Avg MSE (out-of-sample): 1.2904
In [47]:
X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(
        n_estimators=250, bootstrap=True, oob_score=True, n_jobs=-1
    )
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))
Avg MSE (within-sample): 0.1767
Avg MSE (out-of-sample): 1.2843
In [48]:
xcols_ls_pred = [v for v in XCOLS_LS if v != 'ls_year']
i_ls_pred = d[xcols_ls_pred].isnull().sum(1).eq(0) & d['p_prot'].le(5)
X_ls_pred = d.loc[i_ls_pred, xcols_ls_pred]
X_ls_pred['ls_year'] = 2010
In [49]:
model_et_coords_250 = ExtraTreesRegressor(
    n_estimators=250, n_jobs=-1, bootstrap=True, oob_score=True
) 
model_et_coords_250.fit(X_ls.join(coords), y_ls)
y_pred_et_coords_250 = pd.Series(
    model_et_coords_250.predict(X_ls_pred.join(coords)), index=X_ls_pred.index
)
In [50]:
fig, ax = plt.subplots(figsize=(17, 10))
dp_et = shp_c.join(y_pred_et_coords_250.rename('pred1'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp_et.plot(
    'pred1',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Values of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_et_coords_250_pred.png", dpi=200)
Out[50]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [51]:
resid = (y - y_pred_et_coords_250).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)
ax.set_title("Residuals of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_et_coords_250_resid.png", dpi=200)
Out[51]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [52]:
model_rf_coords_100 = RandomForestRegressor(
    n_estimators=100, n_jobs=-1, oob_score=True
) 
model_rf_coords_100.fit(X_ls.join(coords), y_ls)
y_pred_rf_coords_100 = pd.Series(
    model_rf_coords_100.predict(X_ls_pred.join(coords)), index=X_ls_pred.index
)
In [53]:
fig, ax = plt.subplots(figsize=(17, 10))
dp_rf = shp_c.join(y_pred_rf_coords_100.rename('pred2'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp_rf.plot(
    'pred2',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Values of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_rf_coords_100_pred.png", dpi=200)
Out[53]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [54]:
resid = (y - y_pred_rf_coords_100).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)
ax.set_title("Residuals of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_rf_coords_100_resid.png", dpi=200)
Out[54]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [55]:
xcols_fc_pred = [v for v in XCOLS_FC if v != 'fc_year']
i_fc_pred = d[xcols_fc_pred].isnull().sum(1).eq(0) & d['p_prot'].le(5)
X_fc_pred = d.loc[i_fc_pred, xcols_fc_pred]
model_fc_coords_250 = ExtraTreesRegressor(
    n_estimators=250, n_jobs=-1, bootstrap=True, oob_score=True
)
model_fc_coords_250.fit(X_fc.join(coords), y_fc)
y_pred_fc_coords_250 = pd.Series(
    model_fc_coords_250.predict(X_fc_pred.join(coords)), index=X_fc_pred.index
)
In [56]:
fig, ax = plt.subplots(figsize=(17, 10))
dp_fc = shp_c.join(y_pred_fc_coords_250.rename('pred3'), how='inner')

shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp_fc.plot(
    'pred3',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral', 
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Forest Cover Change in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_fc_et_pred.png", dpi=200)
Out[56]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [57]:
yp = pd.read_csv(wf("y_predictions.csv"))
yp = yp.set_index('pid')
In [58]:
dm = d.join(yp, how='inner')
dm
Out[58]:
fips csd_id ha slope travel p_wet coast_2500 river_frontage lake_frontage lake_importance ... hh_inc_med_tract_2012-2016 hh_inc_avg_bg_2012-2016 hh_inc_avg_tract_2012-2016 p_f_change ls_price_per_ha_log ls_year p_e_ls bld_area_per_ha ls_pred fc_pred
pid
F_803174_3027860 25009 16250 1.138 3.269 8.000 0.0 0.000 0.0 0.0 0.000 ... 98264.0 107637.264151 111868.256722 0 NaN NaN 0.0 3124.780316 12.870125 -2.555327
M_238118_942044 25009 27620 3.052 6.315 14.000 10.0 0.000 0.0 0.0 0.000 ... 96512.0 127585.099685 115417.941953 0 9.598621 1991.0 0.0 1471.821756 12.174028 -4.203807
M_240445_943056 25009 27620 2.165 2.538 20.000 48.0 0.000 0.0 0.0 0.000 ... 96512.0 127585.099685 115417.941953 0 NaN NaN 0.0 3155.196305 12.654431 -1.659536
M_245474_951109 25009 77150 2.162 1.784 19.000 21.0 7.000 0.0 0.0 0.000 ... 138947.0 168559.326425 163715.534591 0 11.840504 2002.0 0.0 1055.504163 12.026018 0.174528
F_778757_3113679 25009 29405 7.070 7.645 8.000 4.0 4.133 0.0 0.0 0.000 ... 72132.0 102680.276817 86602.895054 0 11.389646 2002.0 0.0 0.000000 10.537294 -2.453324
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
M_74987_889210 25003 4545 89.442 9.172 41.429 7.0 0.000 0.0 381.2 148.681 ... 72981.0 89094.137931 87422.722723 0 8.566898 2000.0 0.0 0.000000 8.879863 0.089586
F_182440_2984399 25003 53960 90.803 4.407 48.250 36.0 0.000 0.0 1836.0 3245.908 ... 80388.0 88131.578947 94874.911482 0 NaN NaN 0.0 0.000000 10.725508 -1.071956
M_74001_894053 25003 4545 95.959 7.229 152.000 31.0 0.000 0.0 1090.4 3458.145 ... 72981.0 89094.137931 87422.722723 -1 NaN NaN 0.0 23.843517 8.929642 -0.356848
M_67459_886161 25003 51580 109.997 4.771 71.333 34.0 0.000 0.0 1640.0 1523.297 ... 69306.0 90927.272727 84795.126706 0 7.518052 1982.0 0.0 22.909716 9.359960 -0.341437
F_172467_2984485 25003 53960 115.045 2.960 46.857 40.0 0.000 0.0 604.4 390.769 ... 80388.0 88131.578947 94874.911482 -3 NaN NaN 0.0 0.000000 11.404854 -2.046122

177277 rows × 61 columns

In [59]:
dm['p_f_loss_pred_100yr'] = dm['fc_pred'].mul(-10).clip(upper=dm['p_f_2010'], lower=0)
dm['ha_f_loss_pred_100yr'] = dm['p_f_loss_pred_100yr'] * dm['ha'] / 100
dm['tco2_loss_pred_100yr'] = dm['ha_f_loss_pred_100yr'] * 400
dm['acq_cost_pred_usd'] = np.exp(dm['ls_pred']) * dm['ha']
dm['ce_cost_pred_usd'] = dm['acq_cost_pred_usd'] * 0.4
dm['ce_cost_$/tco2'] = dm['ce_cost_pred_usd'] / dm['tco2_loss_pred_100yr']
In [60]:
dm_filtered = dm[~dm['ce_cost_$/tco2'].isin([np.inf, -np.inf])]
In [61]:
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dm0 = shp_c.join(dm_filtered['ce_cost_$/tco2'].rename('avoidedCO2'), how='inner')
dm0.plot(
    'avoidedCO2',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Conservation Easement Cost per Ton of Expected Avoided CO2 Emissions", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/2_ce_cost_usd_per_tco2.png", dpi=200)
Out[61]:
(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)
In [62]:
dm_sorted = dm_filtered.sort_values(by='ce_cost_$/tco2')
dm_sorted['cum_tco2'] = dm_sorted['tco2_loss_pred_100yr'].cumsum()
dm_supply_curve = dm_sorted[dm_sorted['ce_cost_$/tco2'] <= 1000]
dm_sorted
Out[62]:
fips csd_id ha slope travel p_wet coast_2500 river_frontage lake_frontage lake_importance ... bld_area_per_ha ls_pred fc_pred p_f_loss_pred_100yr ha_f_loss_pred_100yr tco2_loss_pred_100yr acq_cost_pred_usd ce_cost_pred_usd ce_cost_$/tco2 cum_tco2
pid
F_727177_2787269 25005 56375 4.738 0.629 22.0 5.0 0.0 0.0 0.0 0.000 ... 0.000000 7.230578 -2.500678 25.006778 1.184821 473.928465 6.543278e+03 2.617311e+03 5.522587e+00 4.739285e+02
F_725608_2786151 25005 56375 2.373 1.725 22.0 6.0 0.0 0.0 0.0 0.000 ... 0.000000 7.399235 -2.711937 27.119369 0.643543 257.417053 3.879222e+03 1.551689e+03 6.027917e+00 7.313455e+02
F_728494_2785591 25005 56375 1.495 1.300 42.0 0.0 0.0 0.0 0.0 0.000 ... 0.000000 7.509981 -2.894286 28.942861 0.432696 173.078310 2.730136e+03 1.092054e+03 6.309598e+00 9.044238e+02
F_727322_2786475 25005 56375 2.578 1.041 22.0 32.0 0.0 0.0 0.0 0.000 ... 0.000000 6.835246 -1.466861 14.668611 0.378157 151.262712 2.397686e+03 9.590746e+02 6.340456e+00 1.055687e+03
F_727190_2782144 25005 56375 1.657 1.097 12.0 2.0 0.0 0.0 0.0 0.000 ... 0.000000 7.896622 -3.852843 38.528434 0.638416 255.366460 4.454324e+03 1.781730e+03 6.977148e+00 1.311053e+03
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
F_689217_3055697 25017 13135 1.198 0.600 5.5 13.0 0.0 0.0 0.0 0.000 ... 16535.058431 14.083366 -0.001182 0.011822 0.000142 0.056649 1.565976e+06 6.263905e+05 1.105746e+07 5.778186e+07
M_120685_935014 25011 47835 1.155 0.000 42.0 48.0 0.0 155.6 0.0 0.000 ... 0.000000 17.268920 -0.028356 0.283560 0.003275 1.310048 3.650722e+07 1.460289e+07 1.114683e+07 5.778186e+07
M_40312_885061 25003 21360 6.050 3.165 86.0 10.0 0.0 0.0 307.6 61.117 ... 0.000000 11.641281 -0.000022 0.000215 0.000013 0.005206 6.878593e+05 2.751437e+05 5.285255e+07 5.778186e+07
M_226411_884436 25021 50250 5.010 4.274 10.0 0.0 0.0 0.0 0.0 0.000 ... 26438.522954 14.671453 -0.000390 0.003901 0.000195 0.078179 1.179150e+07 4.716598e+06 6.033049e+07 5.778186e+07
M_215342_906477 25017 35425 1.550 2.032 12.0 61.0 0.0 0.0 0.0 0.000 ... 3386.451613 13.715193 -0.000046 0.000457 0.000007 0.002831 1.402053e+06 5.608211e+05 1.981098e+08 5.778186e+07

144079 rows × 68 columns

In [63]:
fig, ax = plt.subplots(figsize=(10, 6))
dm_supply_curve.plot('cum_tco2', 'ce_cost_$/tco2', ax=ax)
ax.set_title("Estimated Supply Curve for Avoided CO2 Emissions", fontsize=18)
ax.set_xlabel("Total Avoided CO2 Emissions (tons)")
ax.set_ylabel("Cost per Ton of Avoided CO2 ($/tCO2)")
#plt.savefig("/Users/andywang/desktop/2_tco2_supply_curve.png", dpi=200)
Out[63]:
Text(0, 0.5, 'Cost per Ton of Avoided CO2 ($/tCO2)')
In [64]:
rggi_price = 13
scc_price = 51
uruguay_price = 135
def outcome_variables(dm, carbon_price):
    selected_parcels = dm[dm['ce_cost_$/tco2'] <= carbon_price]
    n = len(selected_parcels)
    ha_pc = selected_parcels['ha'].sum()
    ha_loss_avo = selected_parcels['ha_f_loss_pred_100yr'].sum()
    tco2_avo = selected_parcels['tco2_loss_pred_100yr'].sum()
    budget = selected_parcels['ce_cost_pred_usd'].sum()
    avg_usd_tco2 = budget / tco2_avo
    return n, ha_pc, ha_loss_avo, tco2_avo, budget, avg_usd_tco2
rggi_outcomes = outcome_variables(dm, rggi_price)
scc_outcomes = outcome_variables(dm, scc_price)
uruguay_outcomes = outcome_variables(dm, uruguay_price)
In [65]:
scenario_df = pd.DataFrame({
    'Scenario': [rggi_price, scc_price, uruguay_price],
    'n': [rggi_outcomes[0], scc_outcomes[0], uruguay_outcomes[0]],
    'ha_pc': [rggi_outcomes[1], scc_outcomes[1], uruguay_outcomes[1]],
    'ha_loss_avo': [rggi_outcomes[2], scc_outcomes[2], uruguay_outcomes[2]],
    'tco2_avo': [rggi_outcomes[3], scc_outcomes[3], uruguay_outcomes[3]],
    'budget': [rggi_outcomes[4], scc_outcomes[4], uruguay_outcomes[4]],
    'avg_usd_tco2': [rggi_outcomes[5], scc_outcomes[5], uruguay_outcomes[5]]
})
#scenario_df.to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False)
In [66]:
BUDGET=70000000
dm_sorted = dm.sort_values(by='ce_cost_$/tco2')
dm_sorted['ce_cost_pred_usd_sum'] = dm_sorted['ce_cost_pred_usd'].cumsum()
i_scenario_1 = dm_sorted['ce_cost_pred_usd_sum'].lt(BUDGET)
selected_parcels_scenario_1 = dm_sorted[i_scenario_1]
In [ ]:
 
In [67]:
n_scenario_1 = len(selected_parcels_scenario_1)
ha_pc_scenario_1 = selected_parcels_scenario_1['ha'].sum()
ha_loss_avo_scenario_1 = selected_parcels_scenario_1['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_1 = selected_parcels_scenario_1['tco2_loss_pred_100yr'].sum()
budget_scenario_1 = selected_parcels_scenario_1['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_1 = budget_scenario_1 / tco2_avo_scenario_1
In [68]:
scenario_1_values = {
    'Scenario': 'Scenario 1',
    'n': n_scenario_1,
    'ha_pc': ha_pc_scenario_1,
    'ha_loss_avo': ha_loss_avo_scenario_1,
    'tco2_avo': tco2_avo_scenario_1,
    'budget': budget_scenario_1,
    'avg_usd_tco2': avg_usd_tco2_scenario_1
}
#pd.DataFrame([scenario_1_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
In [69]:
fig1, ax1 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax1, color='none', edgecolor='black')
selected_parcels_scenario_1_gdf = shp_c.join(selected_parcels_scenario_1, on='pid')
selected_parcels_scenario_1_gdf.plot(ax=ax1, markersize='ha', color='blue', alpha=0.7)
ax1.axis('off')
ax1.set_title("Scenario 1", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_1.png", dpi=200)
Out[69]:
Text(0.5, 1.0, 'Scenario 1')
In [70]:
dm2010 = dm[dm['p_f_2010']>0]
dm_sorted2 = dm2010.sort_values(by='fc_pred')
dm_sorted2['ce_cost_pred_usd_sum'] = dm_sorted2['ce_cost_pred_usd'].cumsum()
i_scenario_2 = dm_sorted2['ce_cost_pred_usd_sum'].lt(BUDGET)
selected_parcels_scenario_2 = dm_sorted2[i_scenario_2]
In [71]:
n_scenario_2 = len(selected_parcels_scenario_2)
ha_pc_scenario_2 = selected_parcels_scenario_2['ha'].sum()
ha_loss_avo_scenario_2 = selected_parcels_scenario_2['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_2 = selected_parcels_scenario_2['tco2_loss_pred_100yr'].sum()
budget_scenario_2 = selected_parcels_scenario_2['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_2 = budget_scenario_2 / tco2_avo_scenario_2
scenario_2_values = {
    'Scenario': 'Scenario 2',
    'n': n_scenario_2,
    'ha_pc': ha_pc_scenario_2,
    'ha_loss_avo': ha_loss_avo_scenario_2,
    'tco2_avo': tco2_avo_scenario_2,
    'budget': budget_scenario_2,
    'avg_usd_tco2': avg_usd_tco2_scenario_2
}
#pd.DataFrame([scenario_2_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
In [72]:
fig2, ax2 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax2, color='none', edgecolor='black')
selected_parcels_scenario_2_gdf = shp_c.join(selected_parcels_scenario_2, on='pid')
selected_parcels_scenario_2_gdf.plot(ax=ax2, markersize='ha', color='blue', alpha=0.7)
ax2.axis('off')
ax2.set_title("Scenario 2", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_2.png", dpi=200)
Out[72]:
Text(0.5, 1.0, 'Scenario 2')
In [91]:
dm['offer3'] = dm['p_f_2010']/100 * dm['ha'] * 2680
i_scenario_3 = dm['ce_cost_pred_usd'] < dm['offer3']
dm[i_scenario_3]['offer3'].sum()
selected_parcels_scenario_3 = dm[i_scenario_3]
dm['offer3']
Out[91]:
pid
F_803174_3027860      1402.9264
M_238118_942044       7933.9792
M_240445_943056       3423.2980
M_245474_951109          0.0000
F_778757_3113679     18758.1240
                       ...     
M_74987_889210      232513.4232
F_182440_2984399     82739.6936
M_74001_894053      218594.6020
M_67459_886161      238781.4876
F_172467_2984485     70913.7380
Name: offer3, Length: 177277, dtype: float64
In [74]:
n_scenario_3 = len(selected_parcels_scenario_3)
ha_pc_scenario_3 = selected_parcels_scenario_3['ha'].sum()
ha_loss_avo_scenario_3 = selected_parcels_scenario_3['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_3 = selected_parcels_scenario_3['tco2_loss_pred_100yr'].sum()
budget_scenario_3 = dm[i_scenario_3]['offer3'].sum()
avg_usd_tco2_scenario_3 = budget_scenario_3 / tco2_avo_scenario_3
scenario_3_values = {
    'Scenario': 'Scenario 3',
    'n': n_scenario_3,
    'ha_pc': ha_pc_scenario_3,
    'ha_loss_avo': ha_loss_avo_scenario_3,
    'tco2_avo': tco2_avo_scenario_3,
    'budget': budget_scenario_3,
    'avg_usd_tco2': avg_usd_tco2_scenario_3
}
#pd.DataFrame([scenario_3_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
In [75]:
fig3, ax3 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax3, color='none', edgecolor='black')
selected_parcels_scenario_3_gdf = shp_c.join(selected_parcels_scenario_3, on='pid')
selected_parcels_scenario_3_gdf.plot(ax=ax3, markersize='ha', color='blue', alpha=0.7)
ax3.axis('off')
ax3.set_title("Scenario 3", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_3.png", dpi=200)
Out[75]:
Text(0.5, 1.0, 'Scenario 3')
In [76]:
dm['offer4']= dm['tco2_loss_pred_100yr'] *41
i_scenario_4 = dm['ce_cost_pred_usd'] < dm['offer4']
dm[i_scenario_4]['offer4'].sum()
selected_parcels_scenario_4 = dm[i_scenario_4]
In [104]:
n_scenario_4 = len(selected_parcels_scenario_4)
ha_pc_scenario_4 = selected_parcels_scenario_4['ha'].sum()
ha_loss_avo_scenario_4 = selected_parcels_scenario_4['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_4 = selected_parcels_scenario_4['tco2_loss_pred_100yr'].sum()
budget_scenario_4 = dm[i_scenario_4]['offer4'].sum()
avg_usd_tco2_scenario_4 = budget_scenario_4 / tco2_avo_scenario_4
scenario_4_values = {
    'Scenario': 'Scenario 4',
    'n': n_scenario_4,
    'ha_pc': ha_pc_scenario_4,
    'ha_loss_avo': ha_loss_avo_scenario_4,
    'tco2_avo': tco2_avo_scenario_4,
    'budget': budget_scenario_4,
    'avg_usd_tco2': avg_usd_tco2_scenario_4
}
#pd.DataFrame([scenario_4_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
In [78]:
fig4, ax4 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax4, color='none', edgecolor='black')
selected_parcels_scenario_4_gdf = shp_c.join(selected_parcels_scenario_4, on='pid')
selected_parcels_scenario_4_gdf.plot(ax=ax4, markersize='ha', color='blue', alpha=0.7)
ax4.axis('off')
ax4.set_title("Scenario 4", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_4.png", dpi=200)
Out[78]:
Text(0.5, 1.0, 'Scenario 4')
In [79]:
transaction_cost_per_parcel = 45000
dm['tc'] = transaction_cost_per_parcel
dm['tc+ce_cost_pred_usd'] = (dm['ce_cost_pred_usd']+dm['tc'])
dm['tc+ce_cost_$/tco2'] = (dm['tc+ce_cost_pred_usd'] / dm['tco2_loss_pred_100yr'])
dm['ce_cost_pred_usd_sum_with_transaction'] = (dm['tc+ce_cost_pred_usd']).cumsum()
dm_sorted5 = dm.sort_values(by='tc+ce_cost_$/tco2')
dm_sorted5
Out[79]:
fips csd_id ha slope travel p_wet coast_2500 river_frontage lake_frontage lake_importance ... tco2_loss_pred_100yr acq_cost_pred_usd ce_cost_pred_usd ce_cost_$/tco2 offer3 offer4 tc tc+ce_cost_pred_usd tc+ce_cost_$/tco2 ce_cost_pred_usd_sum_with_transaction
pid
F_354523_3074640 25011 35285 26.969 11.842 125.00 0.0 0.0 0.0 162.0 11.252 ... 8050.989018 2.485530e+05 99421.213240 12.348944 68663.0740 330090.549737 45000 144421.213240 17.938319 4.168910e+10
M_183328_915554 25027 34165 29.479 3.284 64.25 31.0 0.0 360.0 890.8 306.000 ... 9807.050489 3.424065e+05 136962.585638 13.965727 75053.5340 402089.070054 45000 181962.585638 18.554262 5.810776e+10
F_518568_2860349 25027 63345 28.102 4.145 29.00 39.0 0.0 0.0 0.0 0.000 ... 9667.088000 3.491053e+05 139642.130123 14.445108 64769.4896 396350.608000 45000 184642.130123 19.100078 5.676355e+10
F_401884_2970006 25015 52560 15.180 9.100 125.00 0.0 0.0 0.0 0.0 0.000 ... 5004.677643 1.431493e+05 57259.701524 11.441237 40682.4000 205191.783353 45000 102259.701524 20.432825 4.652155e+10
F_506066_2886601 25027 68155 71.482 6.060 54.75 22.0 0.0 0.0 0.0 0.000 ... 3351.532904 1.079113e+05 43164.514806 12.879037 189656.0424 137412.849083 45000 88164.514806 26.305729 5.684580e+10
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
M_139321_889963 25015 72880 2.359 6.500 27.50 148.0 0.0 426.0 0.0 0.000 ... 0.000000 1.289157e+06 515662.866323 inf 2212.7420 0.000000 45000 560662.866323 inf 4.597590e+10
M_244856_915772 25009 37490 2.013 7.229 1.00 0.0 0.0 0.0 0.0 0.000 ... 0.000000 1.589957e+06 635982.937821 inf 0.0000 0.000000 45000 680982.937821 inf 4.897034e+09
M_240297_934723 25009 7420 1.074 0.864 58.00 0.0 0.0 0.0 0.0 0.000 ... 0.000000 3.828205e+05 153128.200192 inf 0.0000 0.000000 45000 198128.200192 inf 4.896196e+09
F_330501_2911704 25015 62745 6.743 3.781 32.50 7.0 0.0 0.0 0.0 0.000 ... 0.000000 3.442470e+05 137698.808514 inf 4517.8100 0.000000 45000 182698.808514 inf 4.597749e+10
M_259920_930818 25009 21850 1.173 1.174 20.00 32.0 2.0 0.0 0.0 0.000 ... 0.000000 4.870326e+05 194813.021562 inf 0.0000 0.000000 45000 239813.021562 inf 4.497034e+09

177277 rows × 73 columns

In [80]:
i_scenario_5 = dm_sorted5['ce_cost_pred_usd_sum_with_transaction'].lt(BUDGET)
selected_parcels_scenario_5 = dm_sorted5[i_scenario_5]
n_scenario_5 = len(selected_parcels_scenario_5)
ha_pc_scenario_5 = selected_parcels_scenario_5['ha'].sum()
ha_loss_avo_scenario_5 = selected_parcels_scenario_5['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_5 = selected_parcels_scenario_5['tco2_loss_pred_100yr'].sum()
budget_scenario_5 = selected_parcels_scenario_5['tc+ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_5 = budget_scenario_5 / tco2_avo_scenario_5

scenario_5_values = {
    'Scenario': 'Scenario 5',
    'n': n_scenario_5,
    'ha_pc': ha_pc_scenario_5,
    'ha_loss_avo': ha_loss_avo_scenario_5,
    'tco2_avo': tco2_avo_scenario_5,
    'budget': budget_scenario_5,
    'avg_usd_tco2': avg_usd_tco2_scenario_5
}
#pd.DataFrame([scenario_5_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
In [81]:
fig5, ax5 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax5, color='none', edgecolor='black')
selected_parcels_scenario_5_gdf = shp_c.join(selected_parcels_scenario_5, on='pid')
selected_parcels_scenario_5_gdf.plot(ax=ax5, markersize='ha', color='blue', alpha=0.7)
ax5.axis('off')
ax5.set_title("Scenario 5", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_5.png", dpi=200)
Out[81]:
Text(0.5, 1.0, 'Scenario 5')
In [82]:
random_noise = np.random.normal(scale=0.25, size=len(dm))
dm['ls_pred_wta'] = dm['ls_pred'] + random_noise + 0.1
dm['wta_ce_cost_pred_usd'] = np.exp(dm['ls_pred_wta']) * dm['ha'] * 0.4
dm_wta = dm[dm['wta_ce_cost_pred_usd'] < dm['ce_cost_pred_usd']]
dm_wta_sorted = dm_wta.sort_values(by='ce_cost_$/tco2')
dm_wta_sorted['ce_cost_pred_usd_sum'] = dm_wta_sorted['ce_cost_pred_usd'].cumsum()
i_scenario_6 = dm_wta_sorted['ce_cost_pred_usd_sum'].lt(BUDGET)
In [83]:
selected_parcels_scenario_6 = dm_wta_sorted[i_scenario_6]
n_scenario_6 = len(selected_parcels_scenario_6)
ha_pc_scenario_6 = selected_parcels_scenario_6['ha'].sum()
ha_loss_avo_scenario_6 = selected_parcels_scenario_6['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_6 = selected_parcels_scenario_6['tco2_loss_pred_100yr'].sum()
budget_scenario_6 = selected_parcels_scenario_6['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_6 = budget_scenario_6 / tco2_avo_scenario_6

scenario_6_values = {
    'Scenario': 'Scenario 6',
    'n': n_scenario_6,
    'ha_pc': ha_pc_scenario_6,
    'ha_loss_avo': ha_loss_avo_scenario_6,
    'tco2_avo': tco2_avo_scenario_6,
    'budget': budget_scenario_6,
    'avg_usd_tco2': avg_usd_tco2_scenario_6
}
#pd.DataFrame([scenario_6_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
In [84]:
n_scenario_6 = len(selected_parcels_scenario_6)
ha_pc_scenario_6 = selected_parcels_scenario_6['ha'].sum()
ha_loss_avo_scenario_6 = selected_parcels_scenario_6['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_6 = selected_parcels_scenario_6['tco2_loss_pred_100yr'].sum()
budget_scenario_6 = selected_parcels_scenario_6['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_6 = budget_scenario_6 / tco2_avo_scenario_6

scenario_6_values = {
    'Scenario': 'Scenario 6',
    'n': n_scenario_6,
    'ha_pc': ha_pc_scenario_6,
    'ha_loss_avo': ha_loss_avo_scenario_6,
    'tco2_avo': tco2_avo_scenario_6,
    'budget': budget_scenario_6,
    'avg_usd_tco2': avg_usd_tco2_scenario_6
}
#pd.DataFrame([scenario_6_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
In [85]:
fig6, ax6 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax6, color='none', edgecolor='black')
selected_parcels_scenario_6_gdf = shp_c.join(selected_parcels_scenario_6, on='pid')
selected_parcels_scenario_6_gdf.plot(ax=ax6, markersize='ha', color='blue', alpha=0.7)
ax6.axis('off')
ax6.set_title("Scenario 6", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_6.png", dpi=200)
Out[85]:
Text(0.5, 1.0, 'Scenario 6')
In [86]:
'''
Carbon offset projects are based on counterfactual thinking
'''
Out[86]:
'\nCarbon offset projects are based on counterfactual thinking\n'
In [87]:
dm.head()
Out[87]:
fips csd_id ha slope travel p_wet coast_2500 river_frontage lake_frontage lake_importance ... ce_cost_pred_usd ce_cost_$/tco2 offer3 offer4 tc tc+ce_cost_pred_usd tc+ce_cost_$/tco2 ce_cost_pred_usd_sum_with_transaction ls_pred_wta wta_ce_cost_pred_usd
pid
F_803174_3027860 25009 16250 1.138 3.269 8.0 0.0 0.000 0.0 0.0 0.0 ... 176858.824298 1520.470734 1402.9264 4769.057131 45000 221858.824298 1907.339658 2.218588e+05 13.264706 262416.453974
M_238118_942044 25009 27620 3.052 6.315 14.0 10.0 0.000 0.0 0.0 0.0 ... 236460.045368 460.755421 7933.9792 21041.232312 45000 281460.045368 548.440400 5.033189e+05 12.033875 205536.951893
M_240445_943056 25009 27620 2.165 2.538 20.0 48.0 0.000 0.0 0.0 0.0 ... 271186.029362 1886.960249 3423.2980 5892.348399 45000 316186.029362 2200.078191 8.195049e+05 12.478441 227423.866295
M_245474_951109 25009 77150 2.162 1.784 19.0 21.0 7.000 0.0 0.0 0.0 ... 144460.400012 inf 0.0000 0.000000 45000 189460.400012 inf 1.008965e+06 12.129188 160160.366651
F_778757_3113679 25009 29405 7.070 7.645 8.0 4.0 4.133 0.0 0.0 0.0 ... 106602.697769 153.650451 18758.1240 28445.804035 45000 151602.697769 218.510632 1.160568e+06 10.710491 126761.272814

5 rows × 75 columns

In [90]:
# Scenario 1: Optimal Allocation
dm_opt = dm.sort_values(by='ce_cost_$/tco2')
dm_opt['ce_cost_pred_usd_sum'] = dm_opt['ce_cost_pred_usd'].cumsum()
i_scenario_opt = dm_opt['ce_cost_pred_usd_sum'].lt(BUDGET*0.8)
selected_parcels_opt = dm_opt[i_scenario_opt]

# Calculating the outcomes for Scenario 1
scenario_opt_results = {
    'n': len(selected_parcels_opt),
    'ha_pc': selected_parcels_opt['ha'].sum(),
    'ha_loss_avo': selected_parcels_opt['ha_f_loss_pred_100yr'].sum(),
    'tco2_avo': selected_parcels_opt['tco2_loss_pred_100yr'].sum(),
    'budget': selected_parcels_opt['ce_cost_pred_usd'].sum(),
    'avg_usd_tco2': selected_parcels_opt['ce_cost_pred_usd'].sum() / selected_parcels_scenario_1['tco2_loss_pred_100yr'].sum()
}

# Scenario 2: Panic Buying
dm_panic = dm[dm['p_f_2010'] > 0].sort_values(by='fc_pred')
dm_panic['ce_cost_pred_usd_sum'] = dm_panic['ce_cost_pred_usd'].cumsum()
i_scenario_panic = dm_panic['ce_cost_pred_usd_sum'].lt(BUDGET*0.2)
selected_parcels_scenario_panic = dm_sorted2[i_scenario_panic]

# Calculating the outcomes for Scenario 2
scenario_panic_results = {
    'n': len(selected_parcels_scenario_panic),
    'ha_pc': selected_parcels_scenario_panic['ha'].sum(),
    'ha_loss_avo': selected_parcels_scenario_panic['ha_f_loss_pred_100yr'].sum(),
    'tco2_avo': selected_parcels_scenario_panic['tco2_loss_pred_100yr'].sum(),
    'budget': selected_parcels_scenario_panic['ce_cost_pred_usd'].sum(),
    'avg_usd_tco2': selected_parcels_scenario_panic['ce_cost_pred_usd'].sum() / selected_parcels_scenario_panic['tco2_loss_pred_100yr'].sum()
}
scenario_panic_results
#scenario_df.to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False)

#pd.DataFrame([scenario_6_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)
Out[90]:
{'n': 9,
 'ha_pc': 53.827,
 'ha_loss_avo': 7.3386,
 'tco2_avo': 2935.4399999999996,
 'budget': 13500891.525943596,
 'avg_usd_tco2': 4599.273541937017}
In [92]:
dm['tax_offer'] = dm['acq_cost_pred_usd'] * 0.5
dm['tax_offer'] = dm['tax_offer'].clip(upper=75000)
In [102]:
dm['flatrate'] = dm['p_f_2010']/100 * dm['ha'] * 1210
i_scenario_tax = dm['tax_offer'] < dm['flatrate']

print(dm[i_scenario_tax]['flatrate'].sum())
selected_parcels_scenario_tax = dm[i_scenario_tax]
68928753.8709
In [103]:
n_scenario_tax = len(selected_parcels_scenario_tax)
ha_pc_scenario_tax = selected_parcels_scenario_tax['ha'].sum()
ha_loss_avo_scenario_tax = selected_parcels_scenario_tax['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_tax = selected_parcels_scenario_tax['tco2_loss_pred_100yr'].sum()
budget_scenario_tax = dm[i_scenario_tax]['flatrate'].sum()
avg_usd_tco2_scenario_tax = budget_scenario_tax / tco2_avo_scenario_tax
scenario_tax_values = {
    'Scenario': 'Scenario Tax Offer',
    'n': n_scenario_tax,
    'ha_pc': ha_pc_scenario_tax,
    'ha_loss_avo': ha_loss_avo_scenario_tax,
    'tco2_avo': tco2_avo_scenario_tax,
    'budget': budget_scenario_tax,
    'avg_usd_tco2': avg_usd_tco2_scenario_tax
}
scenario_tax_values
Out[103]:
{'Scenario': 'Scenario Tax Offer',
 'n': 587,
 'ha_pc': 66944.37700000001,
 'ha_loss_avo': 6057.395807710985,
 'tco2_avo': 2422958.323084394,
 'budget': 68928753.8709,
 'avg_usd_tco2': 28.448179737220826}
In [118]:
fig, ax10 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax10, color='none', edgecolor='black')
selected_parcels_scenario_tax_gdf = shp_c.join(selected_parcels_scenario_tax, on='pid')
selected_parcels_scenario_tax_gdf.plot(ax=ax10, markersize='ha', color='blue', alpha=0.7)
ax10.axis('off')
ax10.set_title("Tax Offer vs Florest Flat Rate", fontsize=18)
Out[118]:
Text(0.5, 1.0, 'Tax Offer vs Florest Flat Rate')
In [112]:
dm['carbonrate']= dm['tco2_loss_pred_100yr'] *13.1
i_scenario_carbonrate = dm['tax_offer'] < dm['carbonrate']
print(dm[i_scenario_carbonrate]['carbonrate'].sum())
69947096.3479845
In [113]:
selected_parcels_scenario_carbonrate = dm[i_scenario_carbonrate]
n_scenario_carbonrate = len(selected_parcels_scenario_carbonrate)
ha_pc_scenario_carbonrate = selected_parcels_scenario_carbonrate['ha'].sum()
ha_loss_avo_scenario_carbonrate = selected_parcels_scenario_carbonrate['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_carbonrate = selected_parcels_scenario_carbonrate['tco2_loss_pred_100yr'].sum()
budget_scenario_carbonrate = dm[i_scenario_carbonrate]['carbonrate'].sum()
avg_usd_tco2_scenario_carbonrate = budget_scenario_carbonrate / tco2_avo_scenario_carbonrate
scenario_carbonrate_values = {
    'Scenario': 'Scenario Tax Offer vs Carbon Flat Rate',
    'n': n_scenario_carbonrate,
    'ha_pc': ha_pc_scenario_carbonrate,
    'ha_loss_avo': ha_loss_avo_scenario_carbonrate,
    'tco2_avo': tco2_avo_scenario_carbonrate,
    'budget': budget_scenario_carbonrate,
    'avg_usd_tco2': avg_usd_tco2_scenario_carbonrate
}
scenario_carbonrate_values
Out[113]:
{'Scenario': 'Scenario Tax Offer vs Carbon Flat Rate',
 'n': 555,
 'ha_pc': 43257.757,
 'ha_loss_avo': 13348.682509157348,
 'tco2_avo': 5339473.003662939,
 'budget': 69947096.3479845,
 'avg_usd_tco2': 13.099999999999998}
In [119]:
fig, ax11 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax11, color='none', edgecolor='black')
selected_parcels_scenario_carbonrate_gdf = shp_c.join(selected_parcels_scenario_carbonrate, on='pid')
selected_parcels_scenario_carbonrate_gdf.plot(ax=ax11, markersize='ha', color='blue', alpha=0.7)
ax11.axis('off')
ax11.set_title("Tax Offer vs Carbon Flat Rate", fontsize=18)
Out[119]:
Text(0.5, 1.0, 'Tax Offer vs Carbon Flat Rate')
In [114]:
dm_sorted_opt = dm.sort_values(by='ce_cost_$/tco2')
dm_sorted_opt['ce_cost_pred_usd_sum'] = dm_sorted_opt['ce_cost_pred_usd'].cumsum()
dm_sorted_pan = dm[dm['p_f_2010'] > 0].sort_values(by='fc_pred')
dm_sorted_pan['ce_cost_pred_usd_sum'] = dm_sorted_pan['ce_cost_pred_usd'].cumsum()
In [122]:
def simulate_allocation(budget_allocation, dm_sorted_opt, dm_sorted_pan, BUDGET):
    """
    Simulates the allocation of budget between the two scenarios and calculates the total CO2 avoided 
    and average cost per tCO2 for the combined allocation.
    """
    # Allocating budget according to the given ratio
    budget_scenario_opt = BUDGET * budget_allocation
    budget_scenario_pan = BUDGET * (1 - budget_allocation)

    # Scenario 1
    i_scenario_opt = dm_sorted_opt['ce_cost_pred_usd_sum'].lt(budget_scenario_opt)
    selected_scenario_opt = dm_sorted_opt[i_scenario_opt]

    # Scenario 2
    i_scenario_pan = dm_sorted_pan['ce_cost_pred_usd_sum'].lt(budget_scenario_pan)
    selected_scenario_pan = dm_sorted_pan[i_scenario_pan]

    # Combined results
    combined_selected = pd.concat([selected_scenario_opt, selected_scenario_pan])
    total_co2_avoided = combined_selected['tco2_loss_pred_100yr'].sum()
    total_budget_used = combined_selected['ce_cost_pred_usd'].sum()
    avg_cost_per_tco2 = total_budget_used / total_co2_avoided if total_co2_avoided > 0 else np.inf

    return total_co2_avoided, avg_cost_per_tco2

allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
results = [simulate_allocation(allocation, dm_sorted_opt, dm_sorted_pan, BUDGET) for allocation in allocations]

# Pairing results with allocations
results_with_allocations = list(zip(allocations, results))
results_with_allocations
Out[122]:
[(0.0, (18688.224000000002, 3661.831554417337)),
 (0.05, (208408.97933571535, 323.76859366971104)),
 (0.1, (349440.7591711177, 199.9041432084736)),
 (0.15000000000000002, (472187.6891205429, 139.24262811426237)),
 (0.2, (590313.2495407679, 117.41693307236335)),
 (0.25, (699898.3108133498, 93.16154606766555)),
 (0.30000000000000004, (803269.6498168104, 85.4234011409276)),
 (0.35000000000000003, (904628.2792642632, 52.2663046392076)),
 (0.4, (1001494.9751070074, 50.588448312530694)),
 (0.45, (1105043.2615579497, 49.196302369869656)),
 (0.5, (1198860.8999801083, 48.22000418811502)),
 (0.55, (1293750.6744031503, 47.432664764532504)),
 (0.6000000000000001, (1380677.2648041507, 46.857144717419835)),
 (0.65, (1472461.0603783657, 46.371398859145835)),
 (0.7000000000000001, (1518150.9516055286, 44.7606240096028)),
 (0.75, (1641723.4929221836, 42.00639230733035)),
 (0.8, (1728889.8029860104, 40.19401911022835)),
 (0.8500000000000001, (1811077.036106914, 37.253067114854645)),
 (0.9, (1893917.2093978468, 35.08576656362387)),
 (0.9500000000000001, (1974863.0185223534, 35.40127655152479)),
 (1.0, (2055343.026564267, 34.05601719336201))]
In [123]:
# Extracting total CO2 avoided and average cost per tCO2 for each allocation
total_co2_avoided = [result[0] for _, result in results_with_allocations]
avg_cost_per_tco2 = [result[1] for _, result in results_with_allocations]

# Plotting the results
fig, ax1 = plt.subplots()

# Bar graph for total CO2 avoided
ax1.bar(allocations, total_co2_avoided, width=0.03, color='b', align='center')
ax1.set_xlabel('Allocation to Scenario 1')
ax1.set_ylabel('Total CO2 Avoided (tCO2)', color='b')
ax1.tick_params('y', colors='b')

# Creating a second y-axis for average cost per tCO2
ax2 = ax1.twinx()
ax2.plot(allocations, avg_cost_per_tco2, 'r-')
ax2.set_ylabel('Average Cost per tCO2 (USD)', color='r')
ax2.tick_params('y', colors='r')

plt.title('Total CO2 Avoided and Average Cost per tCO2 by Allocation')
plt.show()
In [185]:
def simulate_allocation2(budget_allocation, dm_sorted_opt, dm_sorted_pan, BUDGET):
    """
    Simulates the allocation of budget between the two scenarios and calculates the total CO2 avoided 
    and average cost per tCO2 for the combined allocation.
    """
    # Allocating budget according to the given ratio
    budget_scenario_opt = BUDGET * budget_allocation
    budget_scenario_pan = BUDGET * (1 - budget_allocation)

    # Scenario 1
    i_scenario_opt = dm_sorted_opt['ce_cost_pred_usd_sum'].lt(budget_scenario_opt)
    selected_scenario_opt = dm_sorted_opt[i_scenario_opt]

    # Scenario 2
    i_scenario_pan = dm_sorted_pan['ce_cost_pred_usd_sum'].lt(budget_scenario_pan)
    selected_scenario_pan = dm_sorted_pan[i_scenario_pan]

    # Combined results
    combined_selected = pd.concat([selected_scenario_opt, selected_scenario_pan])
    total_co2_avoided = (selected_scenario_opt['tco2_loss_pred_100yr'].sum()) + ((selected_scenario_pan['tco2_loss_pred_100yr']*15).sum())
    total_budget_used = combined_selected['ce_cost_pred_usd'].sum()
    avg_cost_per_tco2 = total_budget_used / total_co2_avoided if total_co2_avoided > 0 else np.inf

    return total_co2_avoided, avg_cost_per_tco2

allocations = np.arange(0, 1.05, 0.001)  # Allocating from 0% to 100% in 5% increments for Scenario 1
results2 = [simulate_allocation2(allocation, dm_sorted_opt, dm_sorted_pan, BUDGET) for allocation in allocations]

# Pairing results with allocations
results_with_allocations2 = list(zip(allocations, results2))
results_with_allocations2
Out[185]:
[(0.0, (280323.36000000004, 244.12210362782247)),
 (0.001, (288713.8439679416, 237.26656636083797)),
 (0.002, (294996.2639833978, 232.4384172849466)),
 (0.003, (300596.9815912618, 228.32142810910585)),
 (0.004, (300596.9815912618, 228.32142810910585)),
 (0.005, (309825.13321610674, 221.88929843743156)),
 (0.006, (316876.0957418652, 217.24652185580013)),
 (0.007, (316876.0957418652, 217.24652185580013)),
 (0.008, (327729.0138672049, 210.51490343465036)),
 (0.009000000000000001, (331758.5987963735, 208.13088225289147)),
 (0.01, (331758.5987963735, 208.13088225289147)),
 (0.011, (341425.6867963736, 202.6468853055562)),
 (0.012, (344588.328306393, 200.9237092670085)),
 (0.013000000000000001, (349668.3489565756, 198.22663018346896)),
 (0.014, (355837.8782911397, 195.06133766776333)),
 (0.015, (359600.0045853751, 193.1934664265863)),
 (0.016, (363904.37898900127, 191.109884838842)),
 (0.017, (367960.4945795689, 189.1992477427866)),
 (0.018000000000000002, (369517.31664440315, 188.47967546312887)),
 (0.019, (373873.61365684855, 186.50146869741778)),
 (0.02, (378433.83365830104, 184.48199597444028)),
 (0.021, (382624.0307919561, 182.67188575369082)),
 (0.022, (386095.6373751653, 181.2055838128344)),
 (0.023, (383237.2373751653, 182.13912526566745)),
 (0.024, (389678.0340913662, 179.45627608916558)),
 (0.025, (392604.841010284, 177.50995823357292)),
 (0.026000000000000002, (396479.51006050914, 175.97124456981916)),
 (0.027, (400538.6202903059, 174.39338756510426)),
 (0.028, (403779.022465399, 173.15779819657237)),
 (0.029, (407463.7217820729, 171.77897739117378)),
 (0.03, (405849.34711978433, 172.0418983990292)),
 (0.031, (409691.6611566987, 170.62518099579214)),
 (0.032, (352433.77304579003, 197.11864800168212)),
 (0.033, (355572.4123710933, 195.56562187792233)),
 (0.034, (355572.4123710933, 195.56562187792233)),
 (0.035, (364465.4886499952, 191.31453042635528)),
 (0.036000000000000004, (366603.39786894445, 190.3250935284417)),
 (0.037, (370574.2855161191, 188.51859195385242)),
 (0.038, (372572.11358663917, 187.62462469407663)),
 (0.039, (373252.25384664215, 178.89315529807138)),
 (0.04, (373252.25384664215, 178.89315529807138)),
 (0.041, (379784.83542128076, 176.1953143073245)),
 (0.042, (380590.11700571456, 175.86930235084364)),
 (0.043000000000000003, (384844.6121527576, 174.1705963298739)),
 (0.044, (387339.25902592135, 173.1924496877539)),
 (0.045, (392393.63683757844, 171.2493666456154)),
 (0.046, (394744.16211920226, 170.3628569292187)),
 (0.047, (397309.9401704534, 169.40770484684194)),
 (0.048, (401084.66194747627, 168.0277480724304)),
 (0.049, (404695.9233357154, 166.73328851817837)),
 (0.05, (404695.9233357154, 166.73328851817837)),
 (0.051000000000000004, (410752.6815293126, 164.61526396354935)),
 (0.052000000000000005, (413017.0259919108, 163.83992612951985)),
 (0.053, (415284.6237039282, 163.07261079617254)),
 (0.054, (419802.61862060556, 161.56956205651335)),
 (0.055, (422404.58482978435, 160.7203730202542)),
 (0.056, (425374.28873072285, 159.76441231572184)),
 (0.057, (428125.21038101305, 158.89117334349015)),
 (0.058, (429877.95166145643, 158.34139604721594)),
 (0.059000000000000004, (433226.24822982715, 157.3037619702776)),
 (0.06, (436960.37935113243, 156.1664978479264)),
 (0.061, (439018.2236581172, 155.54862184847946)),
 (0.062, (439018.2236581172, 155.54862184847946)),
 (0.063, (439018.2236581172, 155.54862184847946)),
 (0.064, (447566.51517171843, 153.04481546025642)),
 (0.065, (450845.5290900913, 152.1100597542526)),
 (0.066, (453454.5242858703, 151.3762939704835)),
 (0.067, (453454.5242858703, 151.3762939704835)),
 (0.068, (457892.50090475427, 150.14751080671147)),
 (0.069, (461388.72695372836, 149.19665907402444)),
 (0.07, (465719.8983784538, 148.03872388287684)),
 (0.07100000000000001, (467501.8613122781, 147.56881998365057)),
 (0.07200000000000001, (470792.55895586975, 146.7105392463034)),
 (0.073, (472007.24949469394, 146.39681257402617)),
 (0.074, (476291.7734618398, 145.303819722673)),
 (0.075, (478230.08559919096, 144.8157929962546)),
 (0.076, (482567.0230449291, 143.73881748944953)),
 (0.077, (482567.0230449291, 143.73881748944953)),
 (0.078, (486488.62774646626, 142.78315777850707)),
 (0.079, (490318.4563113604, 141.86546920231984)),
 (0.08, (493528.78760656324, 141.10748297354905)),
 (0.081, (495739.19889395055, 140.5914506163418)),
 (0.082, (495739.19889395055, 140.5914506163418)),
 (0.083, (495739.19889395055, 140.5914506163418)),
 (0.084, (504869.6820776064, 138.50791285229096)),
 (0.085, (506316.8146928262, 138.1848779785772)),
 (0.08600000000000001, (500013.96579206135, 138.40637202091614)),
 (0.08700000000000001, (502363.68156453595, 137.87884078411696)),
 (0.088, (504409.4666319285, 137.42394059212313)),
 (0.089, (507703.79646929255, 136.69984713311092)),
 (0.09, (511160.47404613445, 135.9507364566432)),
 (0.091, (511160.47404613445, 135.9507364566432)),
 (0.092, (511160.47404613445, 135.9507364566432)),
 (0.093, (519418.33873952867, 134.20253908016448)),
 (0.094, (521588.3468784452, 133.7525741856015)),
 (0.095, (523994.61129577603, 133.2581373110374)),
 (0.096, (525423.6396596634, 132.96666880833)),
 (0.097, (529993.002688166, 132.04589694736475)),
 (0.098, (528735.7809496236, 131.91767334868416)),
 (0.099, (531113.5265943466, 131.4452014856597)),
 (0.1, (532717.5591711176, 131.12887751046026)),
 (0.101, (534843.8527912907, 130.71250884003157)),
 (0.10200000000000001, (536703.5473733388, 129.57211955185755)),
 (0.10300000000000001, (536890.9974453915, 129.536136353385)),
 (0.10400000000000001, (541330.3946336457, 128.69150529805495)),
 (0.105, (541330.3946336457, 128.69150529805495)),
 (0.106, (545760.5855520514, 127.86255053857808)),
 (0.107, (549874.4043234477, 127.10525227415737)),
 (0.108, (552234.5218079343, 126.67647113548112)),
 (0.109, (533160.5680381829, 118.0557292627597)),
 (0.11, (534164.3349295823, 117.88497526442607)),
 (0.111, (537374.247575688, 117.34359258561803)),
 (0.112, (539462.0340519447, 116.99523335522038)),
 (0.113, (539462.0340519447, 116.99523335522038)),
 (0.114, (544480.3200676299, 116.1690487730297)),
 (0.115, (544480.3200676299, 116.1690487730297)),
 (0.116, (544480.3200676299, 116.1690487730297)),
 (0.117, (553638.8633606753, 114.70097611530228)),
 (0.11800000000000001, (553638.8633606753, 114.70097611530228)),
 (0.11900000000000001, (558650.8394436841, 113.91852875624943)),
 (0.12, (558650.8394436841, 113.91852875624943)),
 (0.121, (563834.0327637601, 113.12477221197025)),
 (0.122, (565734.6757205155, 112.83751054717287)),
 (0.123, (568501.0164592904, 112.42319270919165)),
 (0.124, (571343.3848273718, 112.00216559106583)),
 (0.125, (573183.017854233, 111.73205273722888)),
 (0.126, (576490.5459604702, 111.25101503106767)),
 (0.127, (577650.4205263024, 111.08372880180481)),
 (0.128, (581043.2146936544, 110.59859863280779)),
 (0.129, (583318.3695082249, 110.27670640647764)),
 (0.13, (586503.939959044, 109.83104677365)),
 (0.131, (587839.7878977822, 109.64571942507959)),
 (0.132, (587839.7878977822, 109.64571942507959)),
 (0.133, (593874.8641012921, 108.81969396651516)),
 (0.134, (596062.7013105411, 108.52476076996172)),
 (0.135, (598803.8537820401, 108.15869477381833)),
 (0.136, (599097.7907706681, 108.11973601490624)),
 (0.137, (599097.7907706681, 108.11973601490624)),
 (0.138, (599097.7907706681, 108.11973601490624)),
 (0.139, (599097.7907706681, 108.11973601490624)),
 (0.14, (610976.0762745829, 106.57787536706317)),
 (0.14100000000000001, (613282.9938773179, 106.2856983026789)),
 (0.14200000000000002, (614929.1643524455, 106.0786268656114)),
 (0.14300000000000002, (616977.2231716219, 105.82266920156547)),
 (0.14400000000000002, (619693.7303328018, 105.4861570104073)),
 (0.145, (619693.7303328018, 105.4861570104073)),
 (0.146, (624825.6268319181, 104.85868096236175)),
 (0.147, (624825.6268319181, 104.85868096236175)),
 (0.148, (628378.3644604621, 104.43031418155745)),
 (0.149, (632306.9070855033, 103.96241960929284)),
 (0.15, (632735.7691205429, 103.91170849678781)),
 (0.151, (636984.6963416561, 103.41345287018373)),
 (0.152, (639642.8853156375, 103.10525448108808)),
 (0.153, (641992.3274013412, 102.83507379228585)),
 (0.154, (642987.1205005466, 102.72134616239529)),
 (0.155, (646425.7080968248, 102.3314365361484)),
 (0.156, (649398.3208933763, 101.99809797815621)),
 (0.157, (649398.3208933763, 101.99809797815621)),
 (0.158, (653411.1464325943, 101.55317085224753)),
 (0.159, (656330.2862286085, 101.23347619912903)),
 (0.16, (658471.8407853036, 101.00075791454245)),
 (0.161, (660293.4520971515, 100.80412016211739)),
 (0.162, (660293.4520971515, 100.80412016211739)),
 (0.163, (664907.5606641974, 100.31100511905426)),
 (0.164, (668105.1827653403, 99.9733533101342)),
 (0.165, (669928.0652482576, 99.7824125322852)),
 (0.166, (671960.5048018614, 99.57093906815477)),
 (0.167, (674743.9034732745, 99.28369596995266)),
 (0.168, (677395.2729331465, 99.01247044877637)),
 (0.169, (679148.7653936103, 98.8344602705621)),
 (0.17, (682233.0495026774, 98.52380489993566)),
 (0.171, (684413.103045214, 98.30603265122684)),
 (0.17200000000000001, (685379.6518221999, 98.20993049014673)),
 (0.17300000000000001, (688935.1996898229, 97.85886046557499)),
 (0.17400000000000002, (690114.8340929761, 97.74323664249005)),
 (0.17500000000000002, (693083.8552168652, 97.45428383643804)),
 (0.176, (694174.6179645985, 97.34884802025516)),
 (0.177, (694174.6179645985, 97.34884802025516)),
 (0.178, (700369.5121673498, 96.75633131508656)),
 (0.179, (700369.5121673498, 96.75633131508656)),
 (0.18, (703362.4245648722, 96.47403644145813)),
 (0.181, (707144.26660278, 96.12091256249995)),
 (0.182, (709544.2544063726, 95.8989019879502)),
 (0.183, (709544.2544063726, 95.8989019879502)),
 (0.184, (712784.2661111713, 95.60168000061336)),
 (0.185, (712784.2661111713, 95.60168000061336)),
 (0.186, (717875.7683000085, 95.14031459056804)),
 (0.187, (717875.7683000085, 95.14031459056804)),
 (0.188, (717875.7683000085, 95.14031459056804)),
 (0.189, (725056.8748309747, 94.50101248278338)),
 (0.19, (727526.1589984655, 94.28413181243603)),
 (0.191, (729888.3709504185, 94.07835601052633)),
 (0.192, (729888.3709504185, 94.07835601052633)),
 (0.193, (733429.4541148982, 93.77243283171589)),
 (0.194, (736808.9638116169, 93.48332372352617)),
 (0.195, (736808.9638116169, 93.48332372352617)),
 (0.196, (740890.4261470108, 93.13785354576332)),
 (0.197, (742461.3342627496, 93.00594400151391)),
 (0.198, (744529.5018685127, 92.83326807717539)),
 (0.199, (748319.6513457019, 92.51945230889893)),
 (0.2, (750861.3295407679, 92.31101481208233)),
 (0.201, (751999.2416194648, 92.218171297144)),
 (0.202, (754128.0571376926, 92.04527101335131)),
 (0.203, (754128.0571376926, 92.04527101335131)),
 (0.20400000000000001, (758340.8973684066, 91.70610990111467)),
 (0.20500000000000002, (761358.6716209534, 91.46557866707401)),
 (0.20600000000000002, (763355.1834041034, 91.30758103564335)),
 (0.20700000000000002, (766297.167251514, 91.0764761415436)),
 (0.20800000000000002, (766297.167251514, 91.0764761415436)),
 (0.209, (770134.1119280759, 90.77834257497054)),
 (0.21, (769913.3172584383, 90.68768548036982)),
 (0.211, (772201.165660851, 90.51185702382942)),
 (0.212, (774422.7223312644, 90.34229146227237)),
 (0.213, (775325.2001757037, 89.42555530237674)),
 (0.214, (775325.2001757037, 89.42555530237674)),
 (0.215, (775325.2001757037, 89.42555530237674)),
 (0.216, (781235.1653724734, 88.98760113797827)),
 (0.217, (781235.1653724734, 88.98760113797827)),
 (0.218, (781235.1653724734, 88.98760113797827)),
 (0.219, (788354.7364356115, 88.46878250107159)),
 (0.22, (790945.2202629738, 88.28236301729834)),
 (0.221, (793406.1285671224, 88.10648428416984)),
 (0.222, (794739.7702038572, 88.01167560135569)),
 (0.223, (792643.1302038572, 88.01226394510375)),
 (0.224, (797619.2540121548, 87.66054172926613)),
 (0.225, (795255.8800965825, 79.71911136987276)),
 (0.226, (798326.5675279559, 79.53433494785655)),
 (0.227, (801566.5563707317, 79.34109484158655)),
 (0.228, (802244.893170045, 79.30084546764832)),
 (0.229, (804345.2086027751, 79.1767193516618)),
 (0.23, (806955.3557080526, 79.02342336220049)),
 (0.231, (806955.3557080526, 79.02342336220049)),
 (0.232, (812468.4081362969, 78.7033240794976)),
 (0.233, (813858.901621398, 78.62337373871996)),
 (0.234, (817042.1817340539, 78.44154200969506)),
 (0.23500000000000001, (818959.9726436192, 78.33280682401968)),
 (0.23600000000000002, (821275.2867328429, 78.20226493515555)),
 (0.23700000000000002, (822139.6466072472, 78.15374256485434)),
 (0.23800000000000002, (825666.2131558554, 77.95696918680822)),
 (0.23900000000000002, (827332.8946999549, 77.86469687014052)),
 (0.24, (828729.0380940626, 77.78774734111815)),
 (0.241, (831866.8581510512, 77.6159311054035)),
 (0.242, (834083.4293254395, 77.49556047366403)),
 (0.243, (835795.9286086089, 77.40314622842011)),
 (0.244, (837896.0994257838, 77.29040107324549)),
 (0.245, (840472.9301223276, 77.15295835939146)),
 (0.246, (842824.2854294993, 77.02828130442657)),
 (0.247, (842824.2854294993, 77.02828130442657)),
 (0.248, (842824.2854294993, 77.02828130442657)),
 (0.249, (842824.2854294993, 77.02828130442657)),
 (0.25, (851519.0388133498, 76.57328345398466)),
 (0.251, (852783.609417615, 76.50791799487546)),
 (0.252, (854863.867832747, 76.40084492998572)),
 (0.253, (856905.1848165947, 76.29633738439829)),
 (0.254, (856905.1848165947, 76.29633738439829)),
 (0.255, (856905.1848165947, 76.29633738439829)),
 (0.256, (863997.0536635033, 75.93723967245062)),
 (0.257, (863997.0536635033, 75.93723967245062)),
 (0.258, (863997.0536635033, 75.93723967245062)),
 (0.259, (863997.0536635033, 75.93723967245062)),
 (0.26, (871899.8238389111, 75.54452699187877)),
 (0.261, (871899.8238389111, 75.54452699187877)),
 (0.262, (871899.8238389111, 75.54452699187877)),
 (0.263, (878642.0858993162, 75.2151126125004)),
 (0.264, (881578.7291506311, 75.07380129472081)),
 (0.265, (881578.7291506311, 75.07380129472081)),
 (0.266, (881578.7291506311, 75.07380129472081)),
 (0.267, (887468.9493053427, 74.79350247367438)),
 (0.268, (889640.2700941954, 74.69130284775541)),
 (0.269, (892358.3720924594, 74.56435966434904)),
 (0.27, (892358.3720924594, 74.56435966434904)),
 (0.271, (896472.5219775252, 74.37380704562877)),
 (0.272, (897066.0552528873, 74.34647085563816)),
 (0.273, (899306.3332382625, 74.24363565712781)),
 (0.274, (899306.3332382625, 74.24363565712781)),
 (0.275, (899306.3332382625, 74.24363565712781)),
 (0.276, (906987.1030810544, 73.89531502665828)),
 (0.277, (907736.9387018699, 73.86165044762157)),
 (0.278, (907736.9387018699, 73.86165044762157)),
 (0.279, (907736.9387018699, 73.86165044762157)),
 (0.28, (907736.9387018699, 73.86165044762157)),
 (0.281, (917072.9896207985, 73.44732587841968)),
 (0.28200000000000003, (917072.9896207985, 73.44732587841968)),
 (0.28300000000000003, (917072.9896207985, 73.44732587841968)),
 (0.28400000000000003, (924023.064850831, 73.1444631343147)),
 (0.28500000000000003, (924349.6675257566, 73.1303542960508)),
 (0.28600000000000003, (926685.9053248056, 73.02982218715236)),
 (0.28700000000000003, (926685.9053248056, 73.02982218715236)),
 (0.28800000000000003, (932145.8147621285, 72.79708381582655)),
 (0.289, (933578.0553352935, 72.73649013238114)),
 (0.29, (936432.8150581808, 72.61635569334952)),
 (0.291, (938031.0131453939, 72.54947245610957)),
 (0.292, (938031.0131453939, 72.54947245610957)),
 (0.293, (942358.8107445145, 72.3696510354346)),
 (0.294, (944740.823092622, 72.27147678960762)),
 (0.295, (946851.9346267383, 72.1848903530583)),
 (0.296, (948409.5470732339, 72.12132562293677)),
 (0.297, (950659.092952385, 72.03000102693508)),
 (0.298, (951421.1694201616, 71.99918594030302)),
 (0.299, (954890.3778168105, 71.85958421480476)),
 (0.3, (954890.3778168105, 71.85958421480476)),
 (0.301, (958267.082524488, 71.72478344675453)),
 (0.302, (958267.082524488, 71.72478344675453)),
 (0.303, (962160.1305195732, 71.5705451718161)),
 (0.304, (962160.1305195732, 71.5705451718161)),
 (0.305, (962160.1305195732, 71.5705451718161)),
 (0.306, (962160.1305195732, 71.5705451718161)),
 (0.307, (970612.0465799898, 71.24038459230063)),
 (0.308, (973026.1265799899, 71.14715051425618)),
 (0.309, (975159.1229433364, 71.06516227593252)),
 (0.31, (976631.6964190034, 71.00880351556016)),
 (0.311, (980371.0870102594, 70.86650968365086)),
 (0.312, (980921.1244758915, 70.8456955844089)),
 (0.313, (983028.6909914156, 70.7661747394687)),
 (0.314, (986462.3183934498, 70.63739275235534)),
 (0.315, (986462.3183934498, 70.63739275235534)),
 (0.316, (990675.6038521358, 70.48073737547618)),
 (0.317, (992646.962734439, 70.40799103223117)),
 (0.318, (992646.962734439, 70.40799103223117)),
 (0.319, (979591.2951021367, 46.14466033721565)),
 (0.32, (979591.2951021367, 46.14466033721565)),
 (0.321, (983079.6384156563, 46.10161224227952)),
 (0.322, (983079.6384156563, 46.10161224227952)),
 (0.323, (986698.5768213919, 46.05728731830655)),
 (0.324, (988446.9855275659, 46.03602093971247)),
 (0.325, (990899.7425428599, 46.00633991115474)),
 (0.326, (993607.0871734044, 45.97378401361214)),
 (0.327, (993607.0871734044, 45.97378401361214)),
 (0.328, (993607.0871734044, 45.97378401361214)),
 (0.329, (993607.0871734044, 45.97378401361214)),
 (0.33, (993607.0871734044, 45.97378401361214)),
 (0.331, (993607.0871734044, 45.97378401361214)),
 (0.332, (993607.0871734044, 45.97378401361214)),
 (0.333, (1007870.6387927651, 45.807357192879465)),
 (0.334, (1009534.417752848, 45.78828022713917)),
 (0.335, (1009534.417752848, 45.78828022713917)),
 (0.336, (1009534.417752848, 45.78828022713917)),
 (0.337, (1016643.5827728522, 45.70761672267661)),
 (0.338, (1018628.602212385, 45.68544575832863)),
 (0.339, (1020849.3079159656, 45.660835648365115)),
 (0.34, (1020849.3079159656, 45.660835648365115)),
 (0.341, (1023340.4970658841, 45.63339816297698)),
 (0.342, (1025100.5569352238, 45.6141198520125)),
 (0.343, (1025100.5569352238, 45.6141198520125)),
 (0.34400000000000003, (1025100.5569352238, 45.6141198520125)),
 (0.34500000000000003, (1025100.5569352238, 45.6141198520125)),
 (0.34600000000000003, (1034791.6597035888, 45.50932037801288)),
 (0.34700000000000003, (1036008.8061160499, 45.49633947704848)),
 (0.34800000000000003, (1037317.166614573, 45.48250967065748)),
 (0.34900000000000003, (1040263.1352642633, 45.451555117597216)),
 (0.35000000000000003, (1040263.1352642633, 45.451555117597216)),
 (0.35100000000000003, (1040263.1352642633, 45.451555117597216)),
 (0.352, (1040263.1352642633, 45.451555117597216)),
 (0.353, (1048922.6817314504, 45.36168629993604)),
 (0.354, (1050667.817586246, 45.34390434602361)),
 (0.355, (1051763.152508797, 45.33279027867306)),
 (0.356, (1051763.152508797, 45.33279027867306)),
 (0.357, (1051763.152508797, 45.33279027867306)),
 (0.358, (1058104.3306244076, 45.26907841688496)),
 (0.359, (1058104.3306244076, 45.26907841688496)),
 (0.36, (1058104.3306244076, 45.26907841688496)),
 (0.361, (1058104.3306244076, 45.26907841688496)),
 (0.362, (1066227.379428694, 45.18898383462343)),
 (0.363, (1068674.9943617187, 45.16514680669464)),
 (0.364, (1068674.9943617187, 45.16514680669464)),
 (0.365, (1068674.9943617187, 45.16514680669464)),
 (0.366, (1075029.8515041443, 45.103932163778545)),
 (0.367, (1075029.8515041443, 45.103932163778545)),
 (0.368, (1078641.8515041443, 45.069552403520085)),
 (0.369, (1081482.8712782052, 45.042704428904585)),
 (0.37, (1081957.1120458066, 45.03825061401699)),
 (0.371, (1081957.1120458066, 45.03825061401699)),
 (0.372, (1081957.1120458066, 45.03825061401699)),
 (0.373, (1081957.1120458066, 45.03825061401699)),
 (0.374, (1090805.5120458067, 44.956015436480875)),
 (0.375, (1090805.5120458067, 44.956015436480875)),
 (0.376, (1094702.2963103165, 44.92028289544919)),
 (0.377, (1097444.6803458408, 44.89545075971302)),
 (0.378, (1097444.6803458408, 44.89545075971302)),
 (0.379, (1097444.6803458408, 44.89545075971302)),
 (0.38, (1097444.6803458408, 44.89545075971302)),
 (0.381, (1105395.6506130889, 44.824256449122906)),
 (0.382, (1106377.3452057398, 44.815549304503634)),
 (0.383, (1108352.6735358622, 44.79807916314286)),
 (0.384, (1109708.4871474945, 44.78612678709219)),
 (0.385, (1112046.3821364664, 44.765607926589944)),
 (0.386, (1115301.3284623947, 44.73725893842764)),
 (0.387, (1117401.7742044409, 44.71912362342939)),
 (0.388, (1117401.7742044409, 44.71912362342939)),
 (0.389, (1117401.7742044409, 44.71912362342939)),
 (0.39, (1123502.9431471962, 44.66700176528907)),
 (0.391, (1125408.454372053, 44.650876656328414)),
 (0.392, (1126343.550455512, 44.64303380092687)),
 (0.393, (1129433.545948443, 44.61730196562377)),
 (0.394, (1130023.709054299, 44.61243391877287)),
 (0.395, (1132368.7603157607, 44.59316558803536)),
 (0.396, (1133936.20494226, 44.58033596572943)),
 (0.397, (1137129.8311070076, 44.55434673993035)),
 (0.398, (1137129.8311070076, 44.55434673993035)),
 (0.399, (1137129.8311070076, 44.55434673993035)),
 (0.4, (1137129.8311070076, 44.55434673993035)),
 (0.401, (1137129.8311070076, 44.55434673993035)),
 (0.402, (1137129.8311070076, 44.55434673993035)),
 (0.403, (1148952.740620424, 44.459855832893325)),
 (0.404, (1148952.740620424, 44.459855832893325)),
 (0.405, (1152519.6978856905, 44.431740909808866)),
 (0.406, (1154766.542387676, 44.414235755432216)),
 (0.40700000000000003, (1154766.542387676, 44.414235755432216)),
 (0.40800000000000003, (1159123.2668122984, 44.38050932992492)),
 (0.40900000000000003, (1160030.2438560913, 44.373541234952796)),
 (0.41000000000000003, (1160030.2438560913, 44.373541234952796)),
 (0.41100000000000003, (1164241.8147716764, 44.34132972969252)),
 (0.41200000000000003, (1167060.126905465, 44.31999478614795)),
 (0.41300000000000003, (1167766.9162008637, 44.314705772889454)),
 (0.41400000000000003, (1170305.3071389853, 44.295862040254704)),
 (0.41500000000000004, (1171485.0192665055, 44.28714182394274)),
 (0.41600000000000004, (1171485.0192665055, 44.28714182394274)),
 (0.417, (1171485.0192665055, 44.28714182394274)),
 (0.418, (1171485.0192665055, 44.28714182394274)),
 (0.419, (1180727.8294586996, 44.21949054250587)),
 (0.42, (1181179.7544269322, 44.216219737270734)),
 (0.421, (1184174.3461762918, 44.19463699689569)),
 (0.422, (1184174.3461762918, 44.19463699689569)),
 (0.423, (1184174.3461762918, 44.19463699689569)),
 (0.424, (1184174.3461762918, 44.19463699689569)),
 (0.425, (1184174.3461762918, 44.19463699689569)),
 (0.426, (1184174.3461762918, 44.19463699689569)),
 (0.427, (1184174.3461762918, 44.19463699689569)),
 (0.428, (1197475.9583915316, 44.10020823352574)),
 (0.429, (1199978.2431162565, 44.082839449599085)),
 (0.43, (1201860.9138055767, 44.06991114329871)),
 (0.431, (1204107.9480728044, 44.05463689693561)),
 (0.432, (1204640.7040083574, 44.05104580691101)),
 (0.433, (1204640.7040083574, 44.05104580691101)),
 (0.434, (1208984.4167873391, 44.02188789760931)),
 (0.435, (1210757.3620163281, 44.010050049183185)),
 (0.436, (1212837.394016328, 43.99621051371344)),
 (0.437, (1212837.394016328, 43.99621051371344)),
 (0.438, (1212837.394016328, 43.99621051371344)),
 (0.439, (1219184.3146512485, 43.95429088963534)),
 (0.44, (1220854.1659796364, 43.943340724923246)),
 (0.441, (1223334.2889821571, 43.927136134300895)),
 (0.442, (1225346.7760332383, 43.914051895862485)),
 (0.443, (1227629.3001118815, 43.89932285098565)),
 (0.444, (1229092.1788793693, 43.89002324518806)),
 (0.445, (1231498.6897075009, 43.87482387970246)),
 (0.446, (1231498.6897075009, 43.87482387970246)),
 (0.447, (1234850.7228359035, 43.85395961101607)),
 (0.448, (1236819.3197508873, 43.84176778012131)),
 (0.449, (1236819.3197508873, 43.84176778012131)),
 (0.45, (1240678.1175579496, 43.81800698991745)),
 (0.451, (1243148.3521327027, 43.80297590642868)),
 (0.452, (1244569.6761026108, 43.79443376895202)),
 (0.453, (1245252.162696369, 43.790364870139356)),
 (0.454, (1248020.884684029, 43.77391862905835)),
 (0.455, (1248020.884684029, 43.77391862905835)),
 (0.456, (1248020.884684029, 43.77391862905835)),
 (0.457, (1248020.884684029, 43.77391862905835)),
 (0.458, (1248020.884684029, 43.77391862905835)),
 (0.459, (1248020.884684029, 43.77391862905835)),
 (0.46, (1248020.884684029, 43.77391862905835)),
 (0.461, (1262132.262879228, 43.69125862023591)),
 (0.462, (1263396.0615856436, 43.683979878471106)),
 (0.463, (1265431.552214655, 43.67233202459821)),
 (0.464, (1268062.067006189, 43.657384889486856)),
 (0.465, (1268062.067006189, 43.657384889486856)),
 (0.466, (1271401.705185094, 43.6385772395376)),
 (0.467, (1271401.705185094, 43.6385772395376)),
 (0.468, (1274936.1778122378, 43.6188939478953)),
 (0.46900000000000003, (1276501.6635176167, 43.610250956735435)),
 (0.47000000000000003, (1279542.3697777884, 43.59353253522435)),
 (0.47100000000000003, (1280697.379403688, 43.58725132986373)),
 (0.47200000000000003, (1282558.0574897, 43.57718074637724)),
 (0.47300000000000003, (1284245.3855149182, 43.56811820370469)),
 (0.47400000000000003, (1286681.9183029726, 43.55518922816362)),
 (0.47500000000000003, (1289201.8775555291, 43.541917057706)),
 (0.47600000000000003, (1290844.6582439695, 43.533324846756855)),
 (0.47700000000000004, (1292282.6280260945, 43.525841656333846)),
 (0.47800000000000004, (1292282.6280260945, 43.525841656333846)),
 (0.47900000000000004, (1292282.6280260945, 43.525841656333846)),
 (0.48, (1292282.6280260945, 43.525841656333846)),
 (0.481, (1292282.6280260945, 43.525841656333846)),
 (0.482, (1301860.3077055295, 43.47676185873371)),
 (0.483, (1303231.8935995481, 43.469813812312324)),
 (0.484, (1303231.8935995481, 43.469813812312324)),
 (0.485, (1307926.5042226862, 43.44616706746732)),
 (0.486, (1309741.621709506, 43.43707619854991)),
 (0.487, (1311240.0797741916, 43.42966105106292)),
 (0.488, (1312822.0875215845, 43.42185276496366)),
 (0.489, (1314905.2476985177, 43.41162790174518)),
 (0.49, (1316711.9646005728, 43.4028151871498)),
 (0.491, (1318814.693597164, 43.392652833027235)),
 (0.492, (1318814.693597164, 43.392652833027235)),
 (0.493, (1322910.4913792717, 43.37316158440394)),
 (0.494, (1324398.507392078, 43.366133160480665)),
 (0.495, (1324398.507392078, 43.366133160480665)),
 (0.496, (1324398.507392078, 43.366133160480665)),
 (0.497, (1324398.507392078, 43.366133160480665)),
 (0.498, (1331337.7969407768, 43.333659508191154)),
 (0.499, (1334495.7559801082, 43.31904193696803)),
 (0.5, (1334495.7559801082, 43.31904193696803)),
 (0.501, (1334495.7559801082, 43.31904193696803)),
 (0.502, (1334495.7559801082, 43.31904193696803)),
 (0.503, (1334495.7559801082, 43.31904193696803)),
 (0.504, (1334495.7559801082, 43.31904193696803)),
 (0.505, (1334495.7559801082, 43.31904193696803)),
 (0.506, (1347502.5416236625, 43.259704721766575)),
 (0.507, (1348141.2465956733, 43.25683613058865)),
 (0.508, (1350630.4162415706, 43.24574282233502)),
 (0.509, (1350630.4162415706, 43.24574282233502)),
 (0.51, (1350630.4162415706, 43.24574282233502)),
 (0.511, (1356774.4678275536, 43.21859676444767)),
 (0.512, (1358932.271675483, 43.20917065001594)),
 (0.513, (1360181.836540783, 43.20372703317565)),
 (0.514, (1360181.836540783, 43.20372703317565)),
 (0.515, (1364580.1776882664, 43.184666089156984)),
 (0.516, (1364580.1776882664, 43.184666089156984)),
 (0.517, (1367934.678682162, 43.170264850848675)),
 (0.518, (1369226.2668632688, 43.164744871516575)),
 (0.519, (1371225.1358790824, 43.15623425368964)),
 (0.52, (1371225.1358790824, 43.15623425368964)),
 (0.521, (1371225.1358790824, 43.15623425368964)),
 (0.522, (1371225.1358790824, 43.15623425368964)),
 (0.523, (1379684.1710197094, 43.12057509147576)),
 (0.524, (1381474.8082966425, 43.11310198538347)),
 (0.525, (1382983.9538093756, 43.10684997571817)),
 (0.526, (1382983.9538093756, 43.10684997571817)),
 (0.527, (1382983.9538093756, 43.10684997571817)),
 (0.528, (1382983.9538093756, 43.10684997571817)),
 (0.529, (1390178.0820348314, 43.077420885256956)),
 (0.53, (1392637.6176728825, 43.06748290358141)),
 (0.531, (1393580.7680392403, 43.063700851245)),
 (0.532, (1395724.8809742378, 43.05523250733681)),
 (0.533, (1395724.8809742378, 43.05523250733681)),
 (0.534, (1399209.3209597277, 43.04166621512637)),
 (0.535, (1401659.1241802094, 43.032224895914005)),
 (0.536, (1401659.1241802094, 43.032224895914005)),
 (0.537, (1401659.1241802094, 43.032224895914005)),
 (0.538, (1407576.2034714874, 43.00974246840563)),
 (0.539, (1408570.9840311836, 43.00602573388254)),
 (0.54, (1410183.332897933, 43.00005179299807)),
 (0.541, (1413018.9791654502, 42.98958751013153)),
 (0.542, (1414350.474571311, 42.98471479267136)),
 (0.543, (1416539.9722347488, 42.976742655459375)),
 (0.544, (1418769.6295302422, 42.96873708721449)),
 (0.545, (1418769.6295302422, 42.96873708721449)),
 (0.546, (1422103.331780041, 42.956970753208424)),
 (0.547, (1424095.9210660155, 42.94998358416385)),
 (0.548, (1425641.6978273685, 42.94462349600489)),
 (0.549, (1425641.6978273685, 42.94462349600489)),
 (0.55, (1429385.5304031502, 42.931763840189795)),
 (0.551, (1429385.5304031502, 42.931763840189795)),
 (0.552, (1429385.5304031502, 42.931763840189795)),
 (0.553, (1429385.5304031502, 42.931763840189795)),
 (0.554, (1429385.5304031502, 42.931763840189795)),
 (0.555, (1438727.3242142356, 42.900073390182754)),
 (0.556, (1440553.9174341452, 42.893959657063085)),
 (0.557, (1440553.9174341452, 42.893959657063085)),
 (0.558, (1443626.7962431898, 42.8837379280005)),
 (0.559, (1443626.7962431898, 42.8837379280005)),
 (0.56, (1443626.7962431898, 42.8837379280005)),
 (0.561, (1449682.5906128837, 42.863927347089884)),
 (0.562, (1449682.5906128837, 42.863927347089884)),
 (0.5630000000000001, (1453311.2794306476, 42.85217136287191)),
 (0.5640000000000001, (1455330.644484309, 42.84568444836188)),
 (0.5650000000000001, (1455330.644484309, 42.84568444836188)),
 (0.5660000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5670000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5680000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5690000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5700000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5710000000000001, (1467881.957072041, 42.80593128438942)),
 (0.5720000000000001, (1470047.5313234595, 42.79914877141473)),
 (0.5730000000000001, (1471355.1205017003, 42.79507760741269)),
 (0.5740000000000001, (1473560.5992964958, 42.788234306404)),
 (0.5750000000000001, (1474193.7798087061, 42.78627961043839)),
 (0.5760000000000001, (1477525.6732188696, 42.776063627209965)),
 (0.577, (1478571.6090527945, 42.77287199343539)),
 (0.578, (1478571.6090527945, 42.77287199343539)),
 (0.579, (1482588.1332066883, 42.760672404295576)),
 (0.58, (1484449.3384963144, 42.75508867557045)),
 (0.581, (1484449.3384963144, 42.75508867557045)),
 (0.582, (1484449.3384963144, 42.75508867557045)),
 (0.583, (1489495.3457282702, 42.74014324052423)),
 (0.584, (1492086.0520159716, 42.73259212611171)),
 (0.585, (1492086.0520159716, 42.73259212611171)),
 (0.586, (1494417.5546012283, 42.72586929656484)),
 (0.587, (1496489.8826967315, 42.7199550120354)),
 (0.588, (1499091.326135834, 42.712598615773274)),
 (0.589, (1499091.326135834, 42.712598615773274)),
 (0.59, (1502356.9655342377, 42.70342024945375)),
 (0.591, (1504888.1505252013, 42.696377683037234)),
 (0.592, (1504888.1505252013, 42.696377683037234)),
 (0.593, (1504888.1505252013, 42.696377683037234)),
 (0.594, (1510226.1206418383, 42.68190693511459)),
 (0.595, (1511644.406202314, 42.67809243631294)),
 (0.596, (1511644.406202314, 42.67809243631294)),
 (0.597, (1514154.245903524, 42.67141668482116)),
 (0.598, (1516312.1208041508, 42.665750354003485)),
 (0.599, (1516312.1208041508, 42.665750354003485)),
 (0.6, (1516312.1208041508, 42.665750354003485)),
 (0.601, (1521283.1644316472, 42.65278934369559)),
 (0.602, (1521283.1644316472, 42.65278934369559)),
 (0.603, (1521283.1644316472, 42.65278934369559)),
 (0.604, (1528000.55015498, 42.63548013845723)),
 (0.605, (1529954.6982695658, 42.63050807786182)),
 (0.606, (1531489.568651363, 42.626701930636486)),
 (0.607, (1532607.031279209, 42.623948443914564)),
 (0.608, (1535031.7208215182, 42.61801482605009)),
 (0.609, (1535031.7208215182, 42.61801482605009)),
 (0.61, (1537767.604763809, 42.61141040249598)),
 (0.611, (1540508.40605797, 42.604878302856456)),
 (0.612, (1540508.40605797, 42.604878302856456)),
 (0.613, (1544465.101935842, 42.59551300004948)),
 (0.614, (1546129.3594892966, 42.59163146959967)),
 (0.615, (1547217.6763131209, 42.589121661685795)),
 (0.616, (1547217.6763131209, 42.589121661685795)),
 (0.617, (1547217.6763131209, 42.589121661685795)),
 (0.618, (1547217.6763131209, 42.589121661685795)),
 (0.619, (1547217.6763131209, 42.589121661685795)),
 (0.62, (1556131.8537093808, 42.56883223033582)),
 (0.621, (1556131.8537093808, 42.56883223033582)),
 (0.622, (1556131.8537093808, 42.56883223033582)),
 (0.623, (1556131.8537093808, 42.56883223033582)),
 (0.624, (1556131.8537093808, 42.56883223033582)),
 (0.625, (1556131.8537093808, 42.56883223033582)),
 (0.626, (1556131.8537093808, 42.56883223033582)),
 (0.627, (1569155.3380843708, 42.53987870805724)),
 (0.628, (1570973.7676858357, 42.53591159966429)),
 (0.629, (1572268.025468364, 42.5330995075112)),
 (0.63, (1572268.025468364, 42.5330995075112)),
 (0.631, (1572268.025468364, 42.5330995075112)),
 (0.632, (1577819.5551494425, 42.52122845245574)),
 (0.633, (1579207.9675316876, 42.51829998991907)),
 (0.634, (1579207.9675316876, 42.51829998991907)),
 (0.635, (1579207.9675316876, 42.51829998991907)),
 (0.636, (1579207.9675316876, 42.51829998991907)),
 (0.637, (1586054.6946943023, 42.50398151178249)),
 (0.638, (1588310.9521224909, 42.499297114068604)),
 (0.639, (1590811.513966618, 42.4941360326921)),
 (0.64, (1592355.9299020236, 42.49099505858163)),
 (0.641, (1593520.3639331204, 42.48864838581197)),
 (0.642, (1596283.7490394549, 42.48314935315228)),
 (0.643, (1596283.7490394549, 42.48314935315228)),
 (0.644, (1599118.192173848, 42.477589287323795)),
 (0.645, (1599118.192173848, 42.477589287323795)),
 (0.646, (1599118.192173848, 42.477589287323795)),
 (0.647, (1605401.5033365383, 42.46537704140558)),
 (0.648, (1605401.5033365383, 42.46537704140558)),
 (0.649, (1608095.9163783656, 42.46020305128399)),
 (0.65, (1608095.9163783656, 42.46020305128399)),
 (0.651, (1608095.9163783656, 42.46020305128399)),
 (0.652, (1613386.584077212, 42.4502984543921)),
 (0.653, (1613386.584077212, 42.4502984543921)),
 (0.654, (1613386.584077212, 42.4502984543921)),
 (0.655, (1613386.584077212, 42.4502984543921)),
 (0.656, (1613386.584077212, 42.4502984543921)),
 (0.657, (1621424.765186395, 42.43539594106051)),
 (0.658, (1623643.8370462516, 42.43131582900795)),
 (0.659, (1623643.8370462516, 42.43131582900795)),
 (0.66, (1627527.3887331164, 42.42427426801734)),
 (0.661, (1629327.8551135068, 42.42104620273736)),
 (0.662, (1631776.9930122467, 42.4167069548855)),
 (0.663, (1631776.9930122467, 42.4167069548855)),
 (0.664, (1631776.9930122467, 42.4167069548855)),
 (0.665, (1636275.3971363776, 42.40879979372882)),
 (0.666, (1639006.0228085304, 42.40413740995934)),
 (0.667, (1640703.0528782867, 42.40125933598711)),
 (0.668, (1642513.4217077591, 42.3982063220022)),
 (0.669, (1642513.4217077591, 42.3982063220022)),
 (0.67, (1642513.4217077591, 42.3982063220022)),
 (0.671, (1647247.6287427414, 42.39027797395335)),
 (0.672, (1649074.0824638307, 42.38727193226452)),
 (0.673, (1651435.4266256928, 42.38345125399764)),
 (0.674, (1646936.0852463106, 42.46163571202437)),
 (0.675, (1596803.4519518015, 42.45672631658782)),
 (0.676, (1599790.8715948632, 42.451761598726776)),
 (0.677, (1600781.3036055283, 42.450135931485825)),
 (0.678, (1600781.3036055283, 42.450135931485825)),
 (0.679, (1600781.3036055283, 42.450135931485825)),
 (0.68, (1600781.3036055283, 42.450135931485825)),
 (0.681, (1600781.3036055283, 42.450135931485825)),
 (0.682, (1600781.3036055283, 42.450135931485825)),
 (0.683, (1600781.3036055283, 42.450135931485825)),
 (0.684, (1600781.3036055283, 42.450135931485825)),
 (0.685, (1600781.3036055283, 42.450135931485825)),
 (0.686, (1600781.3036055283, 42.450135931485825)),
 (0.687, (1600781.3036055283, 42.450135931485825)),
 (0.6880000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6890000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6900000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6910000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6920000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6930000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6940000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6950000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6960000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6970000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6980000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6990000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7000000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7010000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7020000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7030000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.704, (1600781.3036055283, 42.450135931485825)),
 (0.705, (1600781.3036055283, 42.450135931485825)),
 (0.706, (1597514.3036055283, 42.20825302448651)),
 (0.707, (1597514.3036055283, 42.20825302448651)),
 (0.708, (1597514.3036055283, 42.20825302448651)),
 (0.709, (1597514.3036055283, 42.20825302448651)),
 (0.71, (1597514.3036055283, 42.20825302448651)),
 (0.711, (1597514.3036055283, 42.20825302448651)),
 (0.712, (1597514.3036055283, 42.20825302448651)),
 (0.713, (1661346.4773899356, 42.1171811454383)),
 (0.714, (1658936.367716731, 42.19183232321753)),
 (0.715, (1645166.127716731, 42.29412572907674)),
 (0.716, (1648303.8016055124, 42.28954262265229)),
 (0.717, (1650044.3829549288, 42.28708710394712)),
 (0.718, (1650928.2129560313, 42.285869266164624)),
 (0.719, (1654002.3016276627, 42.28166071596513)),
 (0.72, (1644744.6249621569, 40.71643451371214)),
 (0.721, (1646642.6656229037, 40.715705551533524)),
 (0.722, (1647353.044084335, 40.71543864277567)),
 (0.723, (1650152.881095927, 40.71440298321249)),
 (0.724, (1650621.739738017, 40.714242397562195)),
 (0.725, (1653015.0989600557, 40.71345367576092)),
 (0.726, (1653015.0989600557, 40.71345367576092)),
 (0.727, (1655901.623090563, 40.71253776382341)),
 (0.728, (1655901.623090563, 40.71253776382341)),
 (0.729, (1660441.162711699, 40.71119680175664)),
 (0.73, (1662132.6269994548, 40.7107660569616)),
 (0.731, (1662132.6269994548, 40.7107660569616)),
 (0.732, (1665589.7387862017, 40.709913446061165)),
 (0.733, (1667339.4993075284, 40.70949949576115)),
 (0.734, (1669169.1785453279, 40.70909818635682)),
 (0.735, (1670529.810126635, 40.70880610735869)),
 (0.736, (1670529.810126635, 40.70880610735869)),
 (0.737, (1670529.810126635, 40.70880610735869)),
 (0.738, (1675310.0249780242, 40.707838279338326)),
 (0.739, (1676688.63114592, 40.70756642913803)),
 (0.74, (1678525.395965629, 40.7072159299589)),
 (0.741, (1678525.395965629, 40.7072159299589)),
 (0.742, (1682747.6938063533, 40.70643498813516)),
 (0.743, (1684254.285889074, 40.70617898799609)),
 (0.744, (1686554.0413080982, 40.70580458695191)),
 (0.745, (1686554.0413080982, 40.70580458695191)),
 (0.746, (1689521.9884616295, 40.7054266326454)),
 (0.747, (1691585.3810323426, 40.705185956203955)),
 (0.748, (1691585.3810323426, 40.705185956203955)),
 (0.749, (1694215.876922184, 40.70489601899619)),
 (0.75, (1694215.876922184, 40.70489601899619)),
 (0.751, (1698686.3964483014, 40.704448633109095)),
 (0.752, (1699402.7382996231, 40.70437717842187)),
 (0.753, (1699402.7382996231, 40.70437717842187)),
 (0.754, (1702828.5672207854, 40.70406598825349)),
 (0.755, (1704281.0093980301, 40.70393793106664)),
 (0.756, (1706256.0487849286, 40.70378028452452)),
 (0.757, (1708959.3634961592, 40.70365248965767)),
 (0.758, (1708959.3634961592, 40.70365248965767)),
 (0.759, (1708959.3634961592, 40.70365248965767)),
 (0.76, (1713394.047050823, 40.70347049496915)),
 (0.761, (1715307.131979079, 40.70345814097489)),
 (0.762, (1716594.8555885812, 40.70347088202108)),
 (0.763, (1718227.3410025674, 40.70349885131586)),
 (0.764, (1716133.3410025674, 40.55907762634134)),
 (0.765, (1719537.6446898018, 40.55944496168652)),
 (0.766, (1722365.5342381606, 40.559840068064574)),
 (0.767, (1722365.5342381606, 40.559840068064574)),
 (0.768, (1716897.3742381604, 40.505705368205604)),
 (0.769, (1716897.3742381604, 40.505705368205604)),
 (0.77, (1723572.5322490223, 40.50689139352863)),
 (0.771, (1725194.647061754, 40.50723824293909)),
 (0.772, (1727239.0620245372, 40.50779449400937)),
 (0.773, (1725108.1020245373, 39.67365096254044)),
 (0.774, (1725108.1020245373, 39.67365096254044)),
 (0.775, (1725108.1020245373, 39.67365096254044)),
 (0.776, (1731795.0810302733, 39.67879641073012)),
 (0.777, (1731795.0810302733, 39.67879641073012)),
 (0.778, (1731795.0810302733, 39.67879641073012)),
 (0.779, (1731795.0810302733, 39.67879641073012)),
 (0.78, (1738662.9980616034, 39.684200292089315)),
 (0.781, (1738662.9980616034, 39.684200292089315)),
 (0.782, (1738662.9980616034, 39.684200292089315)),
 (0.783, (1738662.9980616034, 39.684200292089315)),
 (0.784, (1745300.3186967513, 39.68952017967793)),
 (0.785, (1747109.219937373, 39.69099813468368)),
 (0.786, (1748405.3441451762, 39.69208365498057)),
 (0.787, (1750580.270896708, 39.693922843051055)),
 (0.788, (1752296.3701315718, 39.69542084430647)),
 (0.789, (1753661.083809007, 39.69661420480114)),
 (0.79, (1754427.6388631272, 39.69728516127052)),
 (0.791, (1757386.4464879336, 39.69992018752057)),
 (0.792, (1758530.7287586448, 39.70094259381216)),
 (0.793, (1758530.7287586448, 39.70094259381216)),
 (0.794, (1761154.3537125774, 39.70331347222855)),
 (0.795, (1760126.1937125775, 39.623622361900786)),
 (0.796, (1764689.9425301438, 39.62797705876281)),
 (0.797, (1764346.7716519544, 39.25399652144479)),
 (0.798, (1764346.7716519544, 39.25399652144479)),
 (0.799, (1768071.061551408, 39.258466710385235)),
 (0.8, (1769985.9629860106, 39.260780160914685)),
 (0.801, (1770425.5686484114, 39.261316270282684)),
 (0.802, (1773283.933414241, 39.2648652000591)),
 (0.803, (1775261.5170704531, 39.26736515233705)),
 (0.804, (1775261.5170704531, 39.26736515233705)),
 (0.805, (1775261.5170704531, 39.26736515233705)),
 (0.806, (1779991.0873030403, 39.2733453020131)),
 (0.807, (1779991.0873030403, 39.2733453020131)),
 (0.808, (1775803.3873030401, 39.310605732826104)),
 (0.809, (1763024.6624059796, 36.6453410683926)),
 (0.81, (1763024.6624059796, 36.6453410683926)),
 (0.811, (1765714.1884242543, 36.65287082355967)),
 (0.812, (1768180.751185677, 36.659788835880526)),
 (0.8130000000000001, (1769733.8363652243, 36.66416090940172)),
 (0.8140000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8150000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8160000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8170000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8180000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8190000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8200000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8210000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8220000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8230000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8240000000000001, (1787099.4219535603, 36.71296000364977)),
 (0.8250000000000001, (1789950.2297993782, 36.72094460144245)),
 (0.8260000000000001, (1791334.1586912917, 36.724817165135235)),
 (0.8270000000000001, (1792705.260634783, 36.728679885504036)),
 (0.8280000000000001, (1794687.9799656023, 36.734269771511165)),
 (0.8290000000000001, (1796401.6225554068, 36.7391056365538)),
 (0.8300000000000001, (1796401.6225554068, 36.7391056365538)),
 (0.8310000000000001, (1799445.5261900988, 36.74767489458937)),
 (0.8320000000000001, (1800375.3827280044, 36.75028784883137)),
 (0.833, (1800375.3827280044, 36.75028784883137)),
 (0.834, (1800375.3827280044, 36.75028784883137)),
 (0.835, (1800375.3827280044, 36.75028784883137)),
 (0.836, (1800375.3827280044, 36.75028784883137)),
 (0.837, (1800375.3827280044, 36.75028784883137)),
 (0.838, (1810647.0467280042, 36.77900399388079)),
 (0.839, (1812703.6892941808, 36.784724374726736)),
 (0.84, (1814976.1802845176, 36.791057503851505)),
 (0.841, (1814976.1802845176, 36.791057503851505)),
 (0.842, (1814976.1802845176, 36.791057503851505)),
 (0.843, (1814976.1802845176, 36.791057503851505)),
 (0.844, (1820457.991637114, 36.8064294368114)),
 (0.845, (1820457.991637114, 36.8064294368114)),
 (0.846, (1824890.7682174442, 36.81882838226799)),
 (0.847, (1824890.7682174442, 36.81882838226799)),
 (0.848, (1828532.6275552162, 36.82900528819303)),
 (0.849, (1830014.0049747524, 36.83314743439176)),
 (0.85, (1831515.7481069143, 36.83734330212511)),
 (0.851, (1832234.3585115648, 36.83935053679468)),
 (0.852, (1832234.3585115648, 36.83935053679468)),
 (0.853, (1832234.3585115648, 36.83935053679468)),
 (0.854, (1837723.499598943, 36.85463449820539)),
 (0.855, (1837723.499598943, 36.85463449820539)),
 (0.856, (1837723.499598943, 36.85463449820539)),
 (0.857, (1842662.2640486788, 36.86839934015767)),
 (0.858, (1842662.2640486788, 36.86839934015767)),
 (0.859, (1842662.2640486788, 36.86839934015767)),
 (0.86, (1842662.2640486788, 36.86839934015767)),
 (0.861, (1842662.2640486788, 36.86839934015767)),
 (0.862, (1851942.8262612366, 36.894125223002156)),
 (0.863, (1852135.6767601073, 36.894657900393334)),
 (0.864, (1852135.6767601073, 36.894657900393334)),
 (0.865, (1855595.8370782284, 36.904216439234766)),
 (0.866, (1855595.8370782284, 36.904216439234766)),
 (0.867, (1860237.0736111205, 36.91714071901851)),
 (0.868, (1860237.0736111205, 36.91714071901851)),
 (0.869, (1860237.0736111205, 36.91714071901851)),
 (0.87, (1864108.3232707994, 36.92787904377297)),
 (0.871, (1864108.3232707994, 36.92787904377297)),
 (0.872, (1868456.9050696285, 36.93989816318943)),
 (0.873, (1868456.9050696285, 36.93989816318943)),
 (0.874, (1868456.9050696285, 36.93989816318943)),
 (0.875, (1873296.6461996494, 36.953225025194975)),
 (0.876, (1875233.8129750078, 36.958543090014985)),
 (0.877, (1875981.297129798, 36.96059545377948)),
 (0.878, (1875981.297129798, 36.96059545377948)),
 (0.879, (1880047.6159740987, 36.971774175475545)),
 (0.88, (1881562.9510791958, 36.97600346467724)),
 (0.881, (1882826.6764298012, 36.97953703965471)),
 (0.882, (1884706.2278034957, 36.98478987502741)),
 (0.883, (1885736.1269511934, 36.98766655306129)),
 (0.884, (1887776.5309306253, 36.99335701489712)),
 (0.885, (1887776.5309306253, 36.99335701489712)),
 (0.886, (1888923.757798214, 34.66762236519408)),
 (0.887, (1890109.71013542, 34.672420872793325)),
 (0.888, (1890109.71013542, 34.672420872793325)),
 (0.889, (1890109.71013542, 34.672420872793325)),
 (0.89, (1890109.71013542, 34.672420872793325)),
 (0.891, (1890109.71013542, 34.672420872793325)),
 (0.892, (1890109.71013542, 34.672420872793325)),
 (0.893, (1890109.71013542, 34.672420872793325)),
 (0.894, (1890109.71013542, 34.672420872793325)),
 (0.895, (1903501.1583506772, 34.72634320775449)),
 (0.896, (1903501.1583506772, 34.72634320775449)),
 (0.897, (1906667.120098088, 34.73901515274373)),
 (0.898, (1906667.120098088, 34.73901515274373)),
 (0.899, (1909499.4747170806, 34.7503570799983)),
 (0.9, (1911710.3133978467, 34.75920835602772)),
 (0.901, (1911710.3133978467, 34.75920835602772)),
 (0.902, (1911710.3133978467, 34.75920835602772)),
 (0.903, (1911710.3133978467, 34.75920835602772)),
 (0.904, (1918685.6745838332, 34.787058327585925)),
 (0.905, (1919892.7406931673, 34.79190090697813)),
 (0.906, (1921779.8193519516, 34.799503114307655)),
 (0.907, (1922058.9918212378, 34.800630868268556)),
 (0.908, (1922058.9918212378, 34.800630868268556)),
 (0.909, (1926215.3556158002, 34.81744854195251)),
 (0.91, (1928403.22378676, 34.826276161846344)),
 (0.911, (1929527.4684798832, 34.83081403729165)),
 (0.912, (1931133.8798334722, 34.83729547615242)),
 (0.913, (1932816.6012057953, 34.844080641229986)),
 (0.914, (1932816.6012057953, 34.844080641229986)),
 (0.915, (1932816.6012057953, 34.844080641229986)),
 (0.916, (1938453.6191123489, 34.86673528351463)),
 (0.917, (1939338.9517181413, 34.8702929740858)),
 (0.918, (1941606.9948264377, 34.87940643946802)),
 (0.919, (1941803.0662233133, 34.88019517015477)),
 (0.92, (1944889.9160013136, 34.89261192770758)),
 (0.921, (1945871.9564044562, 34.89655438732529)),
 (0.922, (1945871.9564044562, 34.89655438732529)),
 (0.923, (1948372.1057341264, 34.90657756583472)),
 (0.924, (1951479.4615389556, 34.91900330190154)),
 (0.925, (1951479.4615389556, 34.91900330190154)),
 (0.926, (1951479.4615389556, 34.91900330190154)),
 (0.927, (1955338.5714012883, 34.93443831587016)),
 (0.928, (1957912.6680791155, 34.94471323919435)),
 (0.929, (1959411.7011460653, 34.95070896558495)),
 (0.93, (1959411.7011460653, 34.95070896558495)),
 (0.931, (1959411.7011460653, 34.95070896558495)),
 (0.932, (1963852.4948032387, 34.96843387235004)),
 (0.933, (1965400.5284358142, 34.97461571386869)),
 (0.934, (1967023.3363093412, 34.98110335771664)),
 (0.935, (1969058.006592242, 34.989324934657056)),
 (0.936, (1970945.3545428757, 34.996961168495766)),
 (0.937, (1972340.9155269475, 35.002629905345074)),
 (0.9380000000000001, (1973403.669157547, 35.00696090561922)),
 (0.9390000000000001, (1973403.669157547, 35.00696090561922)),
 (0.9400000000000001, (1976253.3243626468, 35.01860190130032)),
 (0.9410000000000001, (1979076.017567865, 35.03010191082997)),
 (0.9420000000000001, (1979832.417567865, 35.03317854622332)),
 (0.9430000000000001, (1979832.417567865, 35.03317854622332)),
 (0.9440000000000001, (1983179.219007039, 35.046783881052946)),
 (0.9450000000000001, (1983179.219007039, 35.046783881052946)),
 (0.9460000000000001, (1986625.9491780247, 35.06076035794786)),
 (0.9470000000000001, (1988100.6967706233, 35.066728732834584)),
 (0.9480000000000001, (1988100.6967706233, 35.066728732834584)),
 (0.9490000000000001, (1992210.0356570226, 35.08335649798272)),
 (0.9500000000000001, (1992656.1225223537, 35.085166517136756)),
 (0.9510000000000001, (1993982.6245453886, 35.08984287267893)),
 (0.9520000000000001, (1991219.6438592963, 35.070661605078186)),
 (0.9530000000000001, (1991219.6438592963, 35.070661605078186)),
 (0.9540000000000001, (1990120.2038592964, 35.00450208191391)),
 (0.9550000000000001, (1995569.2289316377, 35.026942873248075)),
 (0.9560000000000001, (1993006.6534266125, 34.07481011099974)),
 (0.9570000000000001, (1994077.7219043418, 34.07973163457018)),
 (0.9580000000000001, (1996344.4802670411, 34.09013744578911)),
 (0.9590000000000001, (1996344.4802670411, 34.09013744578911)),
 (0.96, (1998247.4800384867, 34.09885676821181)),
 (0.961, (2000844.9055861712, 34.11074420674395)),
 (0.962, (2001760.3819523714, 34.11493664528704)),
 (0.963, (2004013.9959245203, 34.125259836110615)),
 (0.964, (2004625.4759434313, 34.128060315732476)),
 (0.965, (2007378.8800907831, 34.14064947501463)),
 (0.966, (2008654.8347714012, 34.14647274480419)),
 (0.967, (2010770.7938715508, 34.1561266898063)),
 (0.968, (2012577.0010369993, 34.164383459511356)),
 (0.969, (2013847.8634748417, 34.17019617159261)),
 (0.97, (2015630.7282606578, 34.17835952480938)),
 (0.971, (2017313.0928216837, 34.186060916855794)),
 (0.972, (2018813.6293186292, 34.19292657343878)),
 (0.973, (2020706.499873178, 34.20160340128224)),
 (0.974, (2020706.499873178, 34.20160340128224)),
 (0.975, (2023437.2064653006, 34.21410310235987)),
 (0.976, (2024848.8874757807, 34.22055267034169)),
 (0.977, (2026407.0671254068, 34.22766623810515)),
 (0.978, (2028756.7199919492, 34.23839413086611)),
 (0.979, (2028756.7199919492, 34.23839413086611)),
 (0.98, (2028756.7199919492, 34.23839413086611)),
 (0.981, (2033288.7490129452, 34.25905252523253)),
 (0.982, (2034224.0850522218, 34.26330945258805)),
 (0.983, (2036264.8484946967, 34.27260845452638)),
 (0.984, (2036853.5869306035, 34.27529108400344)),
 (0.985, (2039837.4389131628, 34.28887613024939)),
 (0.986, (2033571.818913163, 34.16509082720591)),
 (0.987, (2035863.1859161742, 34.1756699952195)),
 (0.988, (2037827.1135792038, 34.18473400046884)),
 (0.989, (2039600.6913852254, 34.19290503419731)),
 (0.99, (2040269.0443053609, 34.195992266346465)),
 (0.991, (2043380.1520602964, 34.21035401725523)),
 (0.992, (2043380.1520602964, 34.21035401725523)),
 (0.993, (2042852.781002447, 33.997195544076384)),
 (0.994, (2045439.5961530032, 34.00938838201095)),
 (0.995, (2045840.7266019671, 34.01128039611291)),
 (0.996, (2045840.7266019671, 34.01128039611291)),
 (0.997, (2050142.421807609, 34.031545722549545)),
 (0.998, (2050142.421807609, 34.031545722549545)),
 (0.999, (2053765.617620457, 34.04859254038655)),
 ...]
In [186]:
# Extracting total CO2 avoided and average cost per tCO2 for each allocation
total_co2_avoided2 = [result[0] for _, result in results_with_allocations2]
avg_cost_per_tco22 = [result[1] for _, result in results_with_allocations2]

# Plotting the results
fig, ax12 = plt.subplots()

# Bar graph for total CO2 avoided
ax12.bar(allocations, total_co2_avoided2, width=0.03, color='b', align='center')
ax12.set_xlabel('Allocation to Scenario 1')
ax12.set_ylabel('Total CO2 Avoided (tCO2)', color='b')
ax12.tick_params('y', colors='b')

# Creating a second y-axis for average cost per tCO2
ax13 = ax12.twinx()
ax13.plot(allocations, avg_cost_per_tco22, 'r-')
ax13.set_ylabel('Average Cost per tCO2 (USD)', color='r')
ax13.tick_params('y', colors='r')

plt.title('Total CO2 Avoided and Average Cost per tCO2 by Allocation')
plt.show()
In [156]:
dm['tax_offer'] = dm['acq_cost_pred_usd'] * 0.5
dm['tax_offer'] = dm['tax_offer'].clip(upper=75000)
dm['flatrate'] = dm['p_f_2010']/100 * dm['ha'] * 1210
i_scenario_for = dm['tax_offer'] < dm['flatrate']
selected_parcels_scenario_for = dm[i_scenario_for]
dm['carbonrate']= dm['tco2_loss_pred_100yr'] *13.1
i_scenario_carbonrate = dm['tax_offer'] < dm['carbonrate']
In [159]:
def simulate_allocation(budget_allocation, dm, BUDGET):
    """
    Simulates the allocation of budget between the two scenarios and calculates the total CO2 avoided 
    and average cost per tCO2 for the combined allocation.
    """
    # Allocating budget according to the given ratio
    budget_scenario_for = BUDGET * budget_allocation
    budget_scenario_car = BUDGET * (1 - budget_allocation)

    # Scenario 1
    i_scenario_for = dm['flatrate'].sum().lt(budget_scenario_for)
    selected_scenario_for = dm[i_scenario_for]

    # Scenario 2
    i_scenario_car = i_scenario_carbonrate['carbonrate'].sum().lt(budget_scenario_car)
    selected_scenario_car = dm[i_scenario_car]

    # Combined results
    combined_selected = pd.concat([selected_scenario_for, selected_scenario_car])
    total_co2_avoided = combined_selected['tco2_loss_pred_100yr'].sum()
    total_budget_used = combined_selected['ce_cost_pred_usd'].sum()
    avg_cost_per_tco2 = total_budget_used / total_co2_avoided if total_co2_avoided > 0 else np.inf

    return total_co2_avoided, avg_cost_per_tco2

allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
results3 = [simulate_allocation(allocation, dm, BUDGET) for allocation in allocations]

# Pairing results with allocations
results_with_allocations3 = list(zip(allocations, results))
results_with_allocations3
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/anaconda3/lib/python3.11/site-packages/pandas/core/indexes/base.py:3653, in Index.get_loc(self, key)
   3652 try:
-> 3653     return self._engine.get_loc(casted_key)
   3654 except KeyError as err:

File ~/anaconda3/lib/python3.11/site-packages/pandas/_libs/index.pyx:147, in pandas._libs.index.IndexEngine.get_loc()

File ~/anaconda3/lib/python3.11/site-packages/pandas/_libs/index.pyx:176, in pandas._libs.index.IndexEngine.get_loc()

File pandas/_libs/hashtable_class_helper.pxi:7080, in pandas._libs.hashtable.PyObjectHashTable.get_item()

File pandas/_libs/hashtable_class_helper.pxi:7088, in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'flatrate'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
Cell In[159], line 27
     24     return total_co2_avoided, avg_cost_per_tco2
     26 allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
---> 27 results3 = [simulate_allocation(allocation, dm, BUDGET) for allocation in allocations]
     29 # Pairing results with allocations
     30 results_with_allocations3 = list(zip(allocations, results))

Cell In[159], line 27, in <listcomp>(.0)
     24     return total_co2_avoided, avg_cost_per_tco2
     26 allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
---> 27 results3 = [simulate_allocation(allocation, dm, BUDGET) for allocation in allocations]
     29 # Pairing results with allocations
     30 results_with_allocations3 = list(zip(allocations, results))

Cell In[159], line 11, in simulate_allocation(budget_allocation, dm, BUDGET)
      8 budget_scenario_car = BUDGET * (1 - budget_allocation)
     10 # Scenario 1
---> 11 i_scenario_for_sum = i_scenario_for['flatrate'].sum().lt(budget_scenario_for)
     12 selected_scenario_for = dm[i_scenario_for_sum]
     14 # Scenario 2

File ~/anaconda3/lib/python3.11/site-packages/pandas/core/series.py:1007, in Series.__getitem__(self, key)
   1004     return self._values[key]
   1006 elif key_is_scalar:
-> 1007     return self._get_value(key)
   1009 if is_hashable(key):
   1010     # Otherwise index.get_value will raise InvalidIndexError
   1011     try:
   1012         # For labels that don't resolve as scalars like tuples and frozensets

File ~/anaconda3/lib/python3.11/site-packages/pandas/core/series.py:1116, in Series._get_value(self, label, takeable)
   1113     return self._values[label]
   1115 # Similar to Index.get_value, but we do not fall back to positional
-> 1116 loc = self.index.get_loc(label)
   1118 if is_integer(loc):
   1119     return self._values[loc]

File ~/anaconda3/lib/python3.11/site-packages/pandas/core/indexes/base.py:3655, in Index.get_loc(self, key)
   3653     return self._engine.get_loc(casted_key)
   3654 except KeyError as err:
-> 3655     raise KeyError(key) from err
   3656 except TypeError:
   3657     # If we have a listlike key, _check_indexing_error will raise
   3658     #  InvalidIndexError. Otherwise we fall through and re-raise
   3659     #  the TypeError.
   3660     self._check_indexing_error(key)

KeyError: 'flatrate'
In [ ]: